5
$\begingroup$

I'm working my way through a nice e-book on Kalman filters, translating to Mathematica as I go and I've hit an interesting problem. In this section: https://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/03-Gaussians.ipynb#Computational-Properties-of-Gaussians - the author gives a formula to compute a product of Gaussian PDFs that is also Gaussian. Googling around will tell you that the product of two Gaussian PDFs is not itself Gaussian, but is proportional to a Gaussian with mean and standard deviation as in the image below (see proof here: http://www.tina-vision.net/docs/memos/2003-003.pdf for source) -- I have implemented his function and it works fine, but I'm curious if TransformedDistribution can be used to arrive at the same distribution.

enter image description here

To make it concrete

d1 = NormalDistribution[m1, s1];
d2 = NormalDistribution[m2, s2];
(* Works perfectly *)
TransformedDistribution[
 x + y, {x \[Distributed] d1, y \[Distributed] d2}]

(* Not a Normally Distributed *)
TransformedDistribution[
 x*y, {x \[Distributed] d1, y \[Distributed] d2}]

(* This doesn't work, but I'm wondering if there is something like \
this that would produce the normal distribution cited in the question \
*)
TransformedDistribution[
 Normalize[x*y, Total], {x \[Distributed] d1, y \[Distributed] d2}]

Edit to add example multiplication using the proposed function.

MultiplyGaussian[g1_, g2_] := 
 Module[{mean1, var1, mean2, var2, mean, variance},
  {mean1, var1} = g1 /. NormalDistribution[m_, s_] :> {m, s^2};
  {mean2, var2} = g2 /. NormalDistribution[m_, s_] :> {m, s^2};
  mean = (var1*mean2 + var2*mean1) / (var1 + var2);
  variance = (var1 * var2) / (var1 + var2);
  NormalDistribution[mean, Sqrt[variance]]
  ]
z1 = NormalDistribution[3, 0.7];
z2 = NormalDistribution[4.5, 1];
Plot[{Legended[PDF[z1, x], "N(3,0.7)"], 
  Legended[PDF[z2, x], "N(4.5,2)"] , 
  Legended[PDF[MultiplyGaussian[z1, z2], x], "Product"]}, {x, 1, 10}]
$\endgroup$
  • $\begingroup$ I'm not sure your source is right: mathworld.wolfram.com/NormalProductDistribution.html $\endgroup$ – bill s Sep 5 at 23:13
  • $\begingroup$ Edited to add code showing that the proposal seems to work and clarify that the 2nd link is a proof that I at least couldn't find a problem with after a quick look. $\endgroup$ – Dan Sep 5 at 23:34
  • $\begingroup$ I think my terminology is sloppy and that has caused the problem (as pointed out by this MathOverflow answer: math.stackexchange.com/questions/101062/…). The result here is for the product of PDFs (where I sloppily said random variable). Does this clarification help see a way to get there in Mathematica? $\endgroup$ – Dan Sep 5 at 23:41
  • $\begingroup$ Please clean-up the "rv" vs "pdf" confusion in the text and the title. Your title should probably be something like "What is the distribution of a random variable with pdf proportional to the product of two normal pdf's?" $\endgroup$ – JimB Sep 6 at 0:40
4
$\begingroup$
(* Get product of two normal pdf's *)
prod= PDF[NormalDistribution[μ1, σ1], x]*PDF[NormalDistribution[μ2, σ2], x];

(* Normalize so that the pdf integrates to 1 *)
pdf = prod/Integrate[prod, {x, -∞, ∞}, Assumptions -> {σ1 > 0, σ2 > 0}];

(* Construct associated distribution *)
d = ProbabilityDistribution[pdf, {x, -∞, ∞}, Assumptions -> {σ1 > 0, σ2 > 0}];

(* Find mean and variance *)
Mean[d]
(* (μ2 σ1^2+μ1 σ2^2)/(σ1^2+σ2^2) *)

Variance[d]
(* (σ1^2 σ2^2)/(σ1^2+σ2^2) *)

This matches what the article says the mean and variance should be. But is it a normal distribution? If the moment generating function is of the same form as for a normal distribution, then it has a normal distribution. (We could also use the characteristic function to do this for this particular distribution.)

(* The log of the moment generating function will be in the following form *)
logCF = Expectation[Exp[t z], z \[Distributed] NormalDistribution[μ, σ]] /. 
   Power[E, x_] -> x // Expand
(* t μ+(t^2 σ^2)/2 *)

(* So we look to see if the moment generating function of distribution d is of the same form *)
Collect[Expectation[Exp[t z], z \[Distributed] d] /. Power[E, x_] -> x // Expand, t]
(* (t^2 σ1^2 σ2^2)/(2 (σ1^2+σ2^2))+t ((μ2 σ1^2)/(σ1^2+σ2^2)+(μ1 σ2^2)/(σ1^2+σ2^2)) *)

And it is.

$\endgroup$
4
$\begingroup$
d1 = NormalDistribution[m1, s1];
d2 = NormalDistribution[m2, s2];

These distributions require that

assume = And @@
  (DistributionParameterAssumptions /@ {d1, d2})

(* m1 ∈ Reals && s1 > 0 && m2 ∈ Reals && s2 > 0 *)

As pointed out by @JimB, the PDF formed by the product of the normal PDFs is

PDFprod = Assuming[assume, PDF[d1, x]*PDF[d2, x]/
     Integrate[PDF[d1, x]*PDF[d2, x],
      {x, -Infinity, Infinity}] // Simplify]

(* (E^(-((m2 s1^2 + m1 s2^2 - (s1^2 + s2^2) x)^2/(
  2 s1^2 s2^2 (s1^2 + s2^2)))) Sqrt[s1^2 + s2^2])/(Sqrt[2 π] s1 s2) *)

Comparing with the PDF of the expected normal distribution

PDFprod == Assuming[assume, PDF[NormalDistribution[
     (m1*s2^2 + m2*s1^2)/(s1^2 + s2^2),
     Sqrt[s1^2*s2^2/(s1^2 + s2^2)]], x] // Simplify]

(* True *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.