Efficient function for life duration probability?

Question

Sent to me by a colleague asking me for assistance (the function as is breaks):

f = Module[{a, l = # + 1},
    a = Array[x, l]; 
    FullSimplify@
     Probability[
      And @@ Map[If[# == l, Tr[Take[a, #]] > x[1], Tr[Take[a, #]] < x[1]] &, 
                  Range[2, l]], 
      Array[x, l] \[Distributed] ProductDistribution[{NormalDistribution[0, 10], l}]]] &;

For positive integer arguments 1,2,3,... it calculates the probabilities

• x[1] + x[2] > x[1]
• x[1] + x[2] < x[1] && x[1] + x[2] + x[3] > x[1]
• x[1] + x[2] < x[1] && x[1] + x[2] + x[3] < x[1] && x[1] + x[2] + x[3] + x[4] > x[1]
• ...

where x[_] are IID on a normal distribution.

This is in essence calculating the probability that a realization of the RV survives as a maximum when compared to the sum of it and the next 1,2,... realizations.

For example:

f/@{1,2,3}

{1/2, 1/8, 1/16}

All good and well and as expected.

However,

f@4

grinds for considerably more time than the above (as in minutes vs small fractions of a second), and returns the unevaluated (albeit with the obvious simplification):

Probability[ x[2] < 0 && x[2] + x[3] < 0 && x[2] + x[3] + x[4] < 0 && x[2] + x[3] + x[4] + x[5] > 0, {x[1], x[2], x[3], x[4], x[5]} [Distributed] ProductDistribution[{NormalDistribution[0, 10], 5}]]

Is this a bug in Mathematica (we are both on 11.3) probability routines, or is it throwing in the towel because complexity goes crazy?

I plan to ponder this presently, figured the usual suspects might find it interesting.

Ideas for an working/efficient implementation?

Answer 1

ClearAll[proB]
proB[n_Integer, s_] := Module[{v = s^2 SparseArray[{i_, j_} :> Min[i, j], {n, n}], 
    m = ConstantArray[0, n], aa = Array[a, n]},
   NIntegrate[PDF[MultinormalDistribution[m, v], aa], 
    Evaluate[## & @@ Append[Thread[{Most@aa, -∞, 0}], {Last@aa, 0, ∞}]],  
    MinRecursion -> 5, Method -> "AdaptiveQuasiMonteCarlo"]];

proB[#, 10] & /@ Range[10]

{0.5, 0.124995, 0.0625022, 0.038893, 0.0266069, 0.0111523, 0.00551413, 0.00281841, 0.00114446, 0.000464532}

tD1 = MultinormalDistribution[ConstantArray[0, #],
    #2^2 SparseArray[{i_, j_} :> Min[i, j], {#, #}]] &;
tD2 = TransformedDistribution[Accumulate[Array[a, #]], 
    Distributed[Array[a, #], ProductDistribution[{NormalDistribution[0, #2], #}]]] &;
And @@ (tD1[#, s] == tD2[#, s] & /@ Range[2, 50])

True

Answer 2

With the ideas presented by JimB and Kglr, I pondered this a bit, and realized that the desired results must be true for any jump distribution that is unbiased, continuous, and symmetrical.

So, on a lark I rewrote the function as:

f2 = Module[{a, l = # + 1},
    a = Array[x, l]; 
    FullSimplify@
     Probability[
      And @@ Map[
        If[# == l, Tr[Take[a, #]] > x[1], Tr[Take[a, #]] < x[1]] &, 
        Range[2, l]], 
      Array[x, l] \[Distributed] 
       ProductDistribution[{UniformDistribution[{-1, 1}], l}]]] &;

This allows computation of exact values by Mathematica without the failure:

res = f2 /@ Range@7 // AbsoluteTiming

{71.8287, {1/2, 1/8, 1/16, 5/128, 7/256, 21/1024, 33/2048}}

Which allows us to get the closed form:

probFn = FullSimplify[FindSequenceFunction[Last@res, n], n ∈ Integers && 1 <= n]

Pochhammer[1/2, -1 + n]/(2 n!)

Checking with machine precision matches JimB's sim well:

Table[probFn, {n, 1., 10}]

{0.5,0.125,0.0625,0.0390625,0.0273438,0.0205078,0.0161133,0.013092,0.01091,0.00927353}

Answer 3

This is just an extended comment showing that the same thing happens with the Probability function:

pd = ProductDistribution[{NormalDistribution[0, 10], 4}]
Probability[x1 > 0, {x1, x2, x3, x4} \[Distributed] pd]
(* 1/2 *)
Probability[x1 < 0 && x1 + x2 > 0, {x1, x2, x3, x4} \[Distributed] pd]
(* 1/8 *)
FullSimplify[Probability[x1 < 0 && x1 + x2 < 0 && x1 + x2 + x3 > 0, 
  {x1, x2, x3, x4} \[Distributed] pd]]
(* 1/16 *)
Probability[x1 < 0 && x1 + x2 < 0 && x1 + x2 + x3 < 0 && x1 + x2 + x3 + x4 > 0, 
  {x1, x2, x3, x4} \[Distributed] pd]
(* Probability[x1<0&&x1+x2<0&&x1+x2+x3<0&&x1+x2+x3+x4>0,
  {x1,x2,x3,x4}\[Distributed]ProductDistribution[{NormalDistribution[0,10],4}]] *)

Also any positive standard deviation gives the same set of probabilities.

Addition:

I suspect the answer for f@4 is 5/128. Here's why using a brute force approach. We constuct a new set of variables with y1 = x1, y2 = x1 + x2, y3 = x1 + x2 + x3, and y4 = x1 + x2 + x3 + x4.

a = {{1, 0, 0, 0}, {1, 1, 0, 0}, {1, 1, 1, 0}, {1, 1, 1, 1}};
a.{x1, x2, x3, x4}
(* {x1, x1 + x2, x1 + x2 + x3, x1 + x2 + x3 + x4} *)
(Σ = a.DiagonalMatrix[{1, 1, 1, 1}].Transpose[a]);
d = MultinormalDistribution[{0, 0, 0, 0}, Σ];
Integrate[PDF[d, {y1, y2, y3, y4}], {y1, -∞, 0}, {y2, -∞, 0}, {y3, -∞, 0}, {y4, 0, ∞}]

$$\frac{\int _{-\infty }^0\int _{-\infty }^0\left(\text{erf}\left(\frac{\text{y3}}{\sqrt{2}}\right)+1\right) e^{\frac{1}{4} \left(-3 \text{y1}^2+2 \text{y1} \text{y3}-\text{y3}^2\right)} \text{erfc}\left(\frac{\text{y1}+\text{y3}}{2}\right)d\text{y3}d\text{y1}}{8 \sqrt{2} \pi }$$

NIntegrate[PDF[d, {y1, y2, y3, y4}], {y1, -∞, 0}, {y2, -∞, 0}, {y3, -∞, 0}, {y4, 0, ∞}]
(* 0.03906249965228397 *)

5/128 // N
(* 0.0390625 *)

Maybe using simulations is the only possible approach for values larger than 4.

Simulation approach

Here are the results using simulations.

a = {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 1, 0, 0, 0, 0, 0, 0, 0, 
    0}, {1, 1, 1, 0, 0, 0, 0, 0, 0, 0},
   {1, 1, 1, 1, 0, 0, 0, 0, 0, 0}, {1, 1, 1, 1, 1, 0, 0, 0, 0, 0}, {1,
     1, 1, 1, 1, 1, 0, 0, 0, 0},
   {1, 1, 1, 1, 1, 1, 1, 0, 0, 0}, {1, 1, 1, 1, 1, 1, 1, 1, 0, 0}, {1,
     1, 1, 1, 1, 1, 1, 1, 1, 0},
   {1, 1, 1, 1, 1, 1, 1, 1, 1, 1}};
(Σ = 
   a.DiagonalMatrix[{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}].Transpose[a]);
d = MultinormalDistribution[{0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0}, Σ];
nsim = 10000000;
SeedRandom[12345];
x = RandomVariate[d, nsim];
prob = {Length[Select[x, #[[1]] > 0 &]]/nsim // N,
  Length[Select[x, #[[1]] < 0 && #[[2]] > 0 &]]/nsim // N,
  Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] > 0 &]]/nsim // N,
  Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] > 0 &]]/nsim // N,
  Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] > 0 &]]/nsim // N,
  Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] < 0 && #[[6]] > 0 &]]/nsim // N,
  Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] < 0 && #[[6]] < 0 && #[[7]] > 0 &]]/nsim // N,
  Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] < 0 && #[[6]] < 0 && #[[7]] < 0 && #[[8]] > 0 &]]/nsim // N,
  Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] < 0 && #[[6]] < 0 && #[[7]] < 0 && #[[8]] < 0 && #[[9]] > 0 &]]/nsim // N,
  Length[Select[x, #[[1]] < 0 && #[[2]] < 0 && #[[3]] < 0 && #[[4]] < 0 && #[[5]] < 0 && #[[6]] < 0 && #[[7]] < 0 && #[[8]] < 0 && #[[9]] < 0 && #[[10]] > 0 &]]/nsim // N}

{0.499966, 0.125064, 0.0625243, 0.0390055, 0.0273741, 0.0204875, 0.0161639, 0.0130903, 0.0109223, 0.0092713}