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

Question

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

Answer 1

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

Answer 2

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