# How to Plot the PDF of Product of two Normals

I've been looking to plot the distribution of the product of two normally distributed variables, x and y, that are uncorrelated. I've tried symbolically calculating the PDF using the TransformedDistributioncommand on MultiNormalDistribution, but this does not seem to work.

I've also tried generating a simulation multivariate normal data and taking the product, but I've found the SmoothHistogram command created a graph that wasn't smooth enough for my purposes.

Finally, I've tried numerically integrating as follows (here I assume values for the means and variances):

 Plot[
NIntegrate[
Exp[-(x - 1)^2/(2)]/ Sqrt[2 π] Exp[-(y - 1)^2/(2)]/Sqrt[2 π] DiracDelta[x y - z]
, {x, -∞, ∞}
, {y, -∞, ∞}
, MinRecursion -> 3
]
, {z, -2, 2}]


But it just returns zero; I'm guessing my numerical methodology doesn't work with the Dirac Delta function.

Thanks!

Let $Z_1$ and $Z_2$ be independent Gaussian random variables with unit mean and unit standard deviation. Let $W = Z_1 Z_2$. Clearly $$\begin{eqnarray} F_W\left(w\right) &=& \Pr\left(W \leqslant w\right) = \Pr\left(Z_1 Z_2 \leqslant w\right) \\ &=& \mathbb{E}\left(\Pr\left(Z_1 Z_2 \leqslant w \mid Z_2\right) \right) \\ &=& \mathbb{E}\left(\Pr\left(Z_1 Z_2 \leqslant w \mid Z_2, Z_2 > 0\right)\right) + \mathbb{E}\left(\Pr\left(Z_1 Z_2 \leqslant w \mid Z_2, Z_2 < 0\right)\right) \\ &=& \mathbb{E}\left(\Pr\left(Z_1 \leqslant \frac{w}{Z_2}, Z_2 > 0 \mid Z_2\right)\right) + \mathbb{E}\left(\Pr\left(-Z_1 \leqslant \frac{w}{-Z_2}, Z_2 < 0 \mid Z_2\right)\right) \end{eqnarray}$$ Hence $$f_W\left(w\right) = F_W^\prime\left(w\right) = \int_0^\infty f_{Z_2}\left(x\right) f_{Z_1}\left(\frac{w}{x}\right) \frac{\mathcal{d}x}{x} + \int_0^\infty f_{-Z_2}\left(x\right) f_{-Z_1}\left(\frac{w}{x}\right) \frac{\mathcal{d}x}{x}$$ Where $f_Z\left(z\right) = \phi\left(z-1\right)$, and $f_{-Z}\left(z\right) = \phi\left(z+1\right)$ and $\phi(x)$ denotes the pdf of the standard normal distribution.

Putting this into code:

prodPDF[z_Real] :=
NIntegrate[(PDF[NormalDistribution[1, 1], x] PDF[
NormalDistribution[1, 1], z/x] +
PDF[NormalDistribution[-1, 1], x] PDF[NormalDistribution[-1, 1],
z/x])/x, {x, 0, Infinity}]


we can verify the agreement with simulation:

hist = Histogram[
RandomVariate[
TransformedDistribution[
z1 z2, {z1 \[Distributed] NormalDistribution[1, 1],
z2 \[Distributed] NormalDistribution[1, 1]}], 10^6], "Scott",
PDF];

Show[hist, Plot[prodPDF[z], {z, -4, 6}, PlotRange -> 1]]


• Nice! Thanks! I think from here I'll be able to generalize it to arbitrary mean and variance :) – Kevin Aug 28 '15 at 19:42

An exact symbolic solution can be obtained in the case when $\mu=0$, with arbitrary $\sigma$. We then have two independent $N(0,\sigma^2)$ random variables, each with pdf $f(x)$:

The pdf of the product of two Normals can then be derived exactly as:

... where I am using the TransformProduct function from the mathStatica package for Mathematica. Here is a plot of the exact pdf obtained, as $\sigma$ changes:

The general case $N(\mu, \sigma^2)$

If the identical parents are $N(\mu, \sigma^2)$, the iid pdf $f(x)$ is:

The solution is again obtainable with:

TransformProduct[{f, f}, z]


which returns:

While Mathematica does not yield a closed-form solution, ... we do have the functional form. Simply replace Integrate by NIntegrate, and it is then easy to plot the pdf say $g(z)$, for any $\mu$ and $\sigma$ combination. The solution obtained here is the same as that derived by Sasha (by hand) above.

The same approach can be applied for different means $\mu_1$, $\mu_2$ and different variances (just enter the two different parent pdf's).

• In current versions, one can now use MellinConvolve[]: Assuming[σ > 0, Simplify[2 MellinConvolve[PDF[NormalDistribution[0, σ], t], PDF[NormalDistribution[0, σ], t], t, Abs[x]]]] – J. M.'s torpor Mar 19 '18 at 4:46

You can actually do the integral in closed form:

f[z_] = Integrate[Exp[-x^2/2] Exp[-y^2/2] DiracDelta[x y - z]/(4 Pi),
{x, -Infinity, Infinity}, {y, -Infinity, Infinity},
Assumptions -> z ∈ Reals]


and then Plot

Plot[f[z], {z, -2, 2}, PlotRange -> All]


To change the variances:

f[z_] = Integrate[Exp[-x^2/(2 sigX^2)] Exp[-y^2/(2 sigY^2)]
DiracDelta[x y - z]/(4 Pi sigX sigY),
{x, -Infinity, Infinity}, {y, -Infinity, Infinity},
Assumptions -> {z ∈ Reals, sigX > 0, sigY > 0}]


Thanks to @chris for noticing the typo.

Update: this integral is already done at Mathworld and there is a Mathematica file you can download.

• I think this solution only works for the standard normal. It doesn't seem to work when I change the means and variances. – Kevin Aug 28 '15 at 16:20
• I simply integrated the expression the OP presented. – bill s Aug 28 '15 at 20:36
• There is a typo in your expression: remove the z before {x,-Infinity,…. – chris Aug 28 '15 at 21:20