Survival Probability for Random Walks

Question

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.

Answer 1

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

Answer 2

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:

and overlay it over the original plot also:

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

Answer 3

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]

Answer 4

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]

Answer 5

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] &]