3
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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..." $\endgroup$ – JimB Dec 9 '18 at 23:22
2
$\begingroup$
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

$\endgroup$
  • 1
    $\begingroup$ +1, see my answer, interesting problem... $\endgroup$ – ciao Dec 10 '18 at 1:33
2
$\begingroup$

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}

$\endgroup$
  • $\begingroup$ +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). $\endgroup$ – JimB Dec 10 '18 at 5:40
  • $\begingroup$ @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. $\endgroup$ – ciao Dec 10 '18 at 8:42
1
$\begingroup$

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}

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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