# What is the distribution of a random variable with pdf proportional to the product of two normal pdf's

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.

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}]

• I'm not sure your source is right: mathworld.wolfram.com/NormalProductDistribution.html Sep 5 '19 at 23:13
• 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.
– Dan
Sep 5 '19 at 23:34
• 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?
– Dan
Sep 5 '19 at 23:41
• 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?"
– JimB
Sep 6 '19 at 0:40

(* 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.

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 *)