Calculate exact Probability of product of normal random variables

Question

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[20]= 1/2 (1-Erf[2 Sqrt[2]]))
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..

Answer 1

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.