7
$\begingroup$

The Survival Probability for a walker starting at the origin is defined as the probability that the walker stays positive through n steps. Thanks to the Sparre-Andersen Theorem I know this PDF is given by

Plot[Binomial[2 n, n]*2^(-2 n), {n, 0, 100}]

However, I want to validate this empirically.

My attempt to validate this for n=100:

FoldList[
  If[#2 < 0, 0, #1 + #2] &, 
  Prepend[Accumulate[RandomVariate[NormalDistribution[0, 1], 100]], 0]]

I wantFoldList to stop if #2 < 0 evaluates to True, not just substitute in 0.

$\endgroup$
1
  • $\begingroup$ Will, are you attempting to empirically show that the probability for survival when n=100 is Binomial[2 (100), (100)]*2^(-2 (100))? So repeatedly run, and count the times you survive through 100 steps? If so, are you trying to "While" out of the FoldList to save CPU cycles? Not clear to me... $\endgroup$
    – MikeY
    Jan 8, 2019 at 20:36

5 Answers 5

10
$\begingroup$

Something seems odd to me about your code. You are summing twice, once with Accumulate and once with FoldList. If this is really what you want then you could use:

SeedRandom[26]
sum = Prepend[Accumulate[RandomVariate[NormalDistribution[0, 1], 100]], 0];

TakeWhile[sum, NonNegative] // Accumulate
8

{0, 1.10708, 1.23211, 2.28173, 3.30295, 4.05759, 5.26123, 6.62964}

This is equivalent to your FoldList construct up to the appropriate point:

FoldList[If[#2 < 0, 0, #1 + #2] &, sum]
{0, 1.10708, 1.23211, 2.28173, 3.30295, 4.05759, 5.26123, 6.62964, 0, ...

Perhaps you meant to only sum once. In that case TakeWhile[sum, NonNegative] is a direct solution but also sub-optimal as it does not provide early exit behavior, which I suspect is what you're actually after here. It is not clear to me if you need the cumulative sum (walk) itself or only its length; if the latter consider this:

SeedRandom[26]
dist = RandomVariate[NormalDistribution[0, 1], 100];

Module[{i = 0},
  Fold[If[# < 0, Return[i, Fold], i++; # + #2] &, 0, dist]
]
8
$\endgroup$
1
  • $\begingroup$ The ListLinePlot of your bottom answer is much closer to the binomial than the one in mine. I think you have it right :) $\endgroup$
    – Carl Lange
    Jan 9, 2019 at 8:29
7
$\begingroup$

We can do this using an implementation of FoldWhileList.

First, implement FoldWhileList using this great answer.

FoldWhileList[f_, test_, start_, secargs_List] := 
 Module[{tag}, 
  If[# === {}, {start}, Prepend[First@#, start]] &@
   Reap[Fold[If[test[##], Sow[f[##], tag], Return[Null, Fold]] &, 
      start, secargs], _, #2 &][[2]]]

Now we simply run this using the test #2 >= 0 (note that the implementation of NestWhile breaks when test stops evaluating True - our implementation of FoldWhileList also does this, therefore we invert the test you originally used.

FoldWhileList[Plus, #2 >= 0 &, 0, 
 Prepend[Accumulate[RandomVariate[NormalDistribution[0, 1], 100]], 0]]

We can now estimate your PDF:

pdf estimate

and overlay it over the original plot also:

overlaid plots

which doesn't seem like a great match - perhaps there's an issue with your original code, as this answer surmises.

$\endgroup$
5
$\begingroup$

It seems to me that this is a problem to which Catch and Throw can be usefully applied.

SeedRandom[1];
Module[{result = {0}, s},
  Catch[
    Fold[
      If[#2 < 0, Throw[Null], result = {result, s = #1 + #2}; s] &,
      0,
      Accumulate[RandomVariate[NormalDistribution[0, 1], 100]]]];
  result // Flatten]

result

$\endgroup$
3
$\begingroup$

How about the following brute force approach:

n = 100;
SeedRandom[12345];
nsim = 1000000;
Total[Table[If[Min[Accumulate[RandomVariate[NormalDistribution[0, 1], n]]] <= 0, 0, 1],
    {i, nsim}]]/nsim // N
(* 0.056092 *)
Binomial[2 n, n] 2^(-2 n) // N
(* 0.0563485 *)

To get all of the values from 1 to 100 "simultaneously"...

SeedRandom[12345];
nsim = 100000;
n = 100;
z = ConstantArray[0, n];
Do[
 x = Accumulate[RandomVariate[NormalDistribution[0, 1], n]];
 i = Flatten[Position[x, _?NonPositive]];
 If[Length[i] > 0, If[i[[1]] > 1, z[[1 ;; i[[1]] - 1]] = z[[1 ;; i[[1]] - 1]] + 1], z = z + 1],
 {j, nsim}]
z = z/nsim;
ListPlot[{z, Table[Binomial[2 j, j] 2^(-2 j), {j, n}]}, PlotRange -> All, ImageSize -> Large]

Simulation and exact formula

$\endgroup$
2
$\begingroup$

Count number of steps before random walk value either goes negative or over $m$ steps are already taken, for $n$ walks. Then count amount of last successful steps on each integer bin, reverse it, accumulate these values (essentially extend last nonnegative value backwards to every value before it), reverse it again to get the original order, drop the extra value that counted paths that continued over $m$ steps and calculate probabilities:

With[{n = 5000, m = 100}, 
 Table[First@
         NestWhile[# + {1, RandomVariate@NormalDistribution[]} &,
                   {-1, 0}, Last@# >= 0 && First@# <= m &], n] // 
       BinCounts[#, {1, Max@# + 1, 1}] & // Reverse // Accumulate // 
    Reverse // Most@#/n & // 
  ListPlot[{#, Table[Binomial[2 j, j] 2^(-2 j), {j, m}]}, 
    PlotRange -> All] &]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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