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

   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.



3 Answers 3


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[
     z1 z2, {z1 \[Distributed] NormalDistribution[1, 1], 
      z2 \[Distributed] NormalDistribution[1, 1]}], 10^6], "Scott", 

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

enter image description here

  • $\begingroup$ Nice! Thanks! I think from here I'll be able to generalize it to arbitrary mean and variance :) $\endgroup$
    – Arnold
    Aug 28, 2015 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)$:

enter image description here

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

enter image description here

... 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:

enter image description here

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

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

enter image description here

The solution is again obtainable with:

TransformProduct[{f, f}, z]

which returns:

enter image description here

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

  • $\begingroup$ In current versions, one can now use MellinConvolve[]: Assuming[σ > 0, Simplify[2 MellinConvolve[PDF[NormalDistribution[0, σ], t], PDF[NormalDistribution[0, σ], t], t, Abs[x]]]] $\endgroup$ Mar 19, 2018 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]

enter image description here

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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