# Calculate exact Probability of product of normal random variables

Find the value of $P[\Pi_{i=1}^{10}X_i > C]$ for $C=5$, where $X_{10\times 1}$ is a random vector with $10$ dimensional Normal Distribution having location parameter $\mu_{10\times 1} = (1,1,\dots,1)$ and the scatter parameter $\Sigma = I_{10}$, where $I_{10}$ is the $10 \times 10$ identity matrix.

My question is related to this question. In the previously asked question main interest was Cauchy Distribution. But my problem is for normal distribution. How I can handle this problem. I am able to use Monte Carlo simulation. So, I get value of this probability with Monte Carlo Simulation. Now, I am interested to compare that value with exact probability.

Update

I mainly use Probability[] but can't solve the problem with this function

c = 5;
n = 1;
Probability[Product[x[i], {i, 1, n}] > c,
Table[x[i] \[Distributed] NormalDistribution[a, b], {i, 1, n}]]

(* Out= 1/2 (1-Erf[2 Sqrt]))
c = 5;
n = 2;
Probability[Product[x[i], {i, 1, n}] > c,
Table[x[i] \[Distributed] NormalDistribution[1, 1], {i, 1, n}]]


This is taking huge amount of time..

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey Feb 25 '15 at 20:07
• Please show the code that you have developed to attempt to solve this problem. – bbgodfrey Feb 25 '15 at 20:08
• @bbgodfrey I update the question and thanks for welcome. – Argha Feb 25 '15 at 21:08

For this problem you are likely better off performing a zillion simulations and determining the proportion of times that the product is greater than c.

Below is described a brute-force way to obtain the desired probability for 2 and 3 variables.

If c > 0 and with 2 random normal variables either both need to be positive or both need to be negative for their product to potentially be larger than c. So we just need to integrate over the appropriate regions:

(* Set threshold value *)
c = 5;

(* Product of 2 normal random variables *)

(* Approximate result using numerical integration assuming c > 0 *)
NIntegrate[
Exp[-(x1 - 1)^2/2] Exp[-(x2 - 1)^2/2]/(2 \[Pi]), {x1, 0, Infinity},
{x2, c/x1, Infinity}] +
NIntegrate[
Exp[-(x1 - 1)^2/2] Exp[-(x2 - 1)^2/2]/(2 \[Pi]), {x1, -Infinity, 0},
{x2, -Infinity, c/x1}]

(* Check on results using simulations *)
nsim = 100000;
N[Total[Table[
If[Product[
RandomVariate[NormalDistribution[1, 1], 2][[i]], {i, 2}] > c,
1, 0], {k, nsim}]]/nsim]


For the product of 3 normals there are 8 combinations of the sign of the three random variables that need to be examined. But only the combinations with an even number of negative values need to be considered if c is positive. For the case of 3 variables this means either all are positive (1 way for this to happen) or exactly two are negative (3 ways for this to happen).

(* Product of 3 normal variables *)
NIntegrate[
Exp[-(x1 - 1)^2/2] Exp[-(x2 - 1)^2/2] Exp[-(x3 - 1)^2/2]/(2 \[Pi])^(3/2),
{x1, 0, Infinity}, {x2, 0, Infinity}, {x3, c/(x1 x2), Infinity}] +
3 NIntegrate[
Exp[-(x1 - 1)^2/2] Exp[-(x2 - 1)^2/2] Exp[-(x3 - 1)^2/2]/(2 \[Pi])^3/2),
{x1, -Infinity, 0}, {x2, -Infinity, 0}, {x3, c/(x1 x2), Infinity}]

(* Check on results using simulations *)
nsim = 100000;
N[Total[Table[
If[Product[RandomVariate[NormalDistribution[1, 1], 3][[i]], {i, 3}] > c,
1, 0], {k, nsim}]]/nsim]


Again, much above 3 variables and you'll start running into numerical integration issues and simulations are probably your best bet.