# Efficient function for life duration probability?

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?

• This problem seems to be related to the probability that all values from a multivariate normal are positive (as described in the answer I have below). So this paper might be of interest: projecteuclid.org/download/pdf_1/euclid.aoms/1177700391. From the abstract: "The orthant probability, i.e. the probability that all the xi's will be simultaneously positive, is not, in general, given by a closed expression..." – JimB Dec 9 '18 at 23:22

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

• +1, see my answer, interesting problem... – ciao Dec 10 '18 at 1:33

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}

• +1. Very nice! I replaced UniformDistribution[{-1,1}] with ProbabilityDistribution[(3/2) x^2, {x, -1, 1}] and found the exact same set of probabilities (but it took at lot, lot longer to execute). – JimB Dec 10 '18 at 5:40
• @JimB - Thanks! Puzzling that MMA gives up with a Normal - my guess is there's some optimization trickery under the covers when it "knows" it's dealing with a Normal, and in this case it goes haywire. I find the universality of the result interesting. – ciao Dec 10 '18 at 8:42

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.

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}

• +1, see my answer, interesting problem... – ciao Dec 10 '18 at 1:33