(*^
::[ Information =
"This is a Mathematica Notebook file. It contains ASCII text, and can be
transferred by email, ftp, or other text-file transfer utility. It should
be read or edited using a copy of Mathematica or MathReader. If you
received this as email, use your mail application or copy/paste to save
everything from the line containing (*^ down to the line containing ^*)
into a plain text file. On some systems you may have to give the file a
name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
The line below identifies what version of Mathematica created this file,
but it can be opened using any other version as well.";
FrontEndVersion = "NeXT Mathematica Notebook Front End Version 2.2";
NeXTStandardFontEncoding;
fontset = title, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L1, e8, 24, "Times"; ;
fontset = subtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L1, e6, 18, "Times"; ;
fontset = subsubtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, L1, e6, 14, "Times"; ;
fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, L1, a20, 18, "Times"; ;
fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, L1, a15, 14, "Times"; ;
fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, L1, a12, 12, "Times"; ;
fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12;
fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 10, "Times"; ;
fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5, 12, "Courier"; ;
fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ;
fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5, 12, "Courier"; ;
fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ;
fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5, 12, "Courier"; ;
fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, L1, 12, "Courier"; ;
fontset = name, inactive, noPageBreakInGroup, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, B65535, L1, 10, "Times"; ;
fontset = header, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12;
fontset = leftheader, inactive, 12;
fontset = footer, inactive, nohscroll, noKeepOnOnePage, preserveAspect, center, M7, L1, 12;
fontset = leftfooter, inactive, 12;
fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 10, "Times"; ;
fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12;
fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Courier"; ;
fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12;
fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12;
fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12;
fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12;
fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12;
paletteColors = 128; automaticGrouping; currentKernel;
]
:[font = title; inactive; preserveAspect; startGroup]
Generalized Boundary Problem for One-Dimensional Random Walks
:[font = subtitle; inactive; preserveAspect; fontName = "Palatino"]
by
David K. Neal
Department of Mathematics
Western Kentucky University
Bowling Green, KY 42101
nealdk@wkuvx1.wku.edu
;[s]
1:0,0;118,-1;
1:1,16,12,Times,1,18,0,0,0;
:[font = subsubtitle; inactive; preserveAspect]
Mathematica in Education
Vol. 3 No.4
Copyright 1994
TELOS/Springer-Verlag
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Abstract
;[s]
1:0,0;8,-1;
1:1,18,13,Palatino,1,18,0,0,0;
:[font = text; inactive; preserveAspect; endGroup]
The boundary problem for non-symmetric one-dimensional random walks
is described and simulated graphically. We estimate the probability of a
particle reaching the top boundary before the lower boundary and estimate
the average time to reach either boundary. Confidence intervals for these
estimates are also given. A matrix solution to these problems is then given
followed by an application to the Gambler's Ruin problem.
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Introduction
:[font = text; inactive; pageBreakBelow; preserveAspect]
A famous problem in probability theory is the ÒDrunkard's WalkÓ which can
be described as follows. Suppose a particle begins at height m > 0 on the y-axis.
On each unit time interval, the particle moves up 1 step with probability p or
down 1 step with probability q = 1-p. (i) What is the probability that the particle
reaches either height n > m or height 0? (ii) What is the probability P[m] that
the particle reaches height n before height 0? (iii) What is the average number
of steps T[m] needed to reach one of these two boundaries?
:[font = text; inactive; pageBreakBelow; preserveAspect; fontName = "Palatino"]
Each problem can be solved by one standard technique involving difference
equations [Chu75, p. 243; Fel57, p. 313]. Most importantly, these solutions yield
closed-form formulas. Namely, (i) the particle will reach a boundary with
probability 1; (ii) for p = .5, P[m] = m/n, and for p .5, P[m] = (1-(q/p)m)/(1-(q/p)n),
and (iii) for p = .5, T[m] = m(n-m), and for p .5, T[m] =
n/(p-q)*(1-(q/p)m)/(1-(q/p)n) - m/(p-q). These formulas can easily be
generalized for any lower boundary less than m.
;[s]
1:0,0;514,-1;
1:1,11,8,Times,0,12,0,0,0;
:[font = text; inactive; pageBreakBelow; preserveAspect; fontName = "Palatino"]
We can generalize the problem further by assuming that the particle moves
up height h with probability p or down height k with probability q = 1-p, where
h and k are any positive real numbers. Unfortunately, the above difference
equation technique will no longer work. This problem has been considered in
connection with sequential sampling when h and k are integers [Fel57, p. 330];
however, there is not a known closed-form solution. In this article, we shall
describe a probabilistic method of approximating the solutions for given values
of p, m, n, h, and k. Then we shall give an algebraic method for finding the
exact solution for any given set of initial values provided h and k are integers.
The techniques will be implemented with brief Mathematica programs.
;[s]
1:0,0;777,-1;
1:1,11,8,Times,0,12,0,0,0;
:[font = text; inactive; preserveAspect; fontName = "Palatino"; endGroup]
The material in this article is one of a dozen or so projects that the author
uses as supplemental material for a beginning calculus-based course in
probability and statistics. For all such projects, students are to read the problem
and solution and then run the program. Students then work a few exercises
which involve slight modifications of the program. These students generally
have had no prior experience with Mathematica or Macintosh machines. By
the end of the semester however, they seem to have a good working knowledge
of both. Moreover, they gain a much greater understanding of probabilistic
concepts and have fun doing it.
;[s]
1:0,0;644,-1;
1:1,11,8,Times,0,12,0,0,0;
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Probability Estimation Technique
:[font = text; inactive; preserveAspect]
To approximate the probability of reaching n first, we can apply the Law
of Large Numbers to a sample of random walks and merely count the
proportion of particles which reach or exceed the top boundary first. We can
also average the number of steps needed to reach or exceed a boundary for this
sample in order to approximate the true average time. We can then establish
bounds for our approximations by constructing confidence intervals.
Mathematica can easily generate random walks which fit the parameters of the
problem and can also display graphs which illustrate the problem.
:[font = text; inactive; preserveAspect]
We shall demonstrate with the following values : p = .40, h = 3.2, k = 2.4,
m = 6, and n = 10. We first enter all the variables:
:[font = input; preserveAspect]
p = .40;
q = 1-p;
StepsUp = h = 3.2;
StepsDown = k = 2.4;
StartingHeight = m = 6;
UpperBoundary = n = 10;
:[font = text; inactive; preserveAspect]
We shall make our estimations based on, say 300, trials. Therefore,
we must construct a loop of 300 random walks which begin at starting height
m and go up h steps with probability p or down k steps with probability q. Each
walk stops when reaching the upper boundary of 10 or the lower boundary of 0.
(If desired, one can always change the lower boundary to any value below m.)
The height of the ith particle on its jth step is denoted by x[i,j], for i = 1,...,300.
For the program, we first set the number of trials and set all the initial heights
equal to m. Within the loop, t[i] counts the number of steps that the ith walk
needed to reach a boundary and u[i] determines whether or not the walk ended
at or above height n.
:[font = input; preserveAspect]
trials = 300;
Do[x[i,0] = m,
{i, 1, trials}];
:[font = input; preserveAspect]
Do[j=0;
While[0 < x[i,j] < n,
x[i,j+1] = If[Random[] < p,
x[i,j] + h,
x[i,j] - k];
j++];
t[i] = j;
u[i] = If[x[i,t[i]] >= n,
1, 0],
{i, 1, trials}]
:[font = text; inactive; preserveAspect]
In order to see graphs of the generated random walks, where the x-axis
denotes time and the y-axis denotes height, we can use the command below.
Here we display only the the first 2 graphs; but the index limits can be changed,
for instance to {i,10,20}, to display any consecutive sequence of graphs.
;[s]
3:0,0;246,1;255,2;304,-1;
3:1,11,8,Times,0,12,0,0,0;1,10,8,Courier,1,12,0,0,0;1,11,8,Times,0,12,0,0,0;
:[font = input; Cclosed; preserveAspect; startGroup]
Do[ListPlot[Table[{j, x[i,j]},
{j, 0, t[i]}],
PlotJoined -> True,
AxesOrigin -> Automatic],
{i, 1, 2}]
:[font = postscript; PostScript; formatAsPostScript; output; inactive; Cclosed; preserveAspect; pictureLeft = 34; pictureWidth = 281; pictureHeight = 173; startGroup]
%!
%%Creator: Mathematica
%%AspectRatio: .61803
MathPictureStart
%% Graphics
/Helvetica findfont 6 scalefont setfont
% Scaling calculations
0.02381 0.190476 0.051503 0.091969 [
[(1)] .21429 .0515 0 2 Msboxa
[(2)] .40476 .0515 0 2 Msboxa
[(3)] .59524 .0515 0 2 Msboxa
[(4)] .78571 .0515 0 2 Msboxa
[(5)] .97619 .0515 0 2 Msboxa
[(1)] .01131 .14347 1 0 Msboxa
[(2)] .01131 .23544 1 0 Msboxa
[(3)] .01131 .32741 1 0 Msboxa
[(4)] .01131 .41938 1 0 Msboxa
[(5)] .01131 .51135 1 0 Msboxa
[(6)] .01131 .60332 1 0 Msboxa
[ -0.001 -0.001 0 0 ]
[ 1.001 .61903 0 0 ]
] MathScale
% Start of Graphics
1 setlinecap
1 setlinejoin
newpath
[ ] 0 setdash
0 g
p
p
.002 w
.21429 .0515 m
.21429 .05775 L
s
P
[(1)] .21429 .0515 0 2 Mshowa
p
.002 w
.40476 .0515 m
.40476 .05775 L
s
P
[(2)] .40476 .0515 0 2 Mshowa
p
.002 w
.59524 .0515 m
.59524 .05775 L
s
P
[(3)] .59524 .0515 0 2 Mshowa
p
.002 w
.78571 .0515 m
.78571 .05775 L
s
P
[(4)] .78571 .0515 0 2 Mshowa
p
.002 w
.97619 .0515 m
.97619 .05775 L
s
P
[(5)] .97619 .0515 0 2 Mshowa
p
.001 w
.0619 .0515 m
.0619 .05525 L
s
P
p
.001 w
.1 .0515 m
.1 .05525 L
s
P
p
.001 w
.1381 .0515 m
.1381 .05525 L
s
P
p
.001 w
.17619 .0515 m
.17619 .05525 L
s
P
p
.001 w
.25238 .0515 m
.25238 .05525 L
s
P
p
.001 w
.29048 .0515 m
.29048 .05525 L
s
P
p
.001 w
.32857 .0515 m
.32857 .05525 L
s
P
p
.001 w
.36667 .0515 m
.36667 .05525 L
s
P
p
.001 w
.44286 .0515 m
.44286 .05525 L
s
P
p
.001 w
.48095 .0515 m
.48095 .05525 L
s
P
p
.001 w
.51905 .0515 m
.51905 .05525 L
s
P
p
.001 w
.55714 .0515 m
.55714 .05525 L
s
P
p
.001 w
.63333 .0515 m
.63333 .05525 L
s
P
p
.001 w
.67143 .0515 m
.67143 .05525 L
s
P
p
.001 w
.70952 .0515 m
.70952 .05525 L
s
P
p
.001 w
.74762 .0515 m
.74762 .05525 L
s
P
p
.001 w
.82381 .0515 m
.82381 .05525 L
s
P
p
.001 w
.8619 .0515 m
.8619 .05525 L
s
P
p
.001 w
.9 .0515 m
.9 .05525 L
s
P
p
.001 w
.9381 .0515 m
.9381 .05525 L
s
P
p
.002 w
0 .0515 m
1 .0515 L
s
P
p
.002 w
.02381 .14347 m
.03006 .14347 L
s
P
[(1)] .01131 .14347 1 0 Mshowa
p
.002 w
.02381 .23544 m
.03006 .23544 L
s
P
[(2)] .01131 .23544 1 0 Mshowa
p
.002 w
.02381 .32741 m
.03006 .32741 L
s
P
[(3)] .01131 .32741 1 0 Mshowa
p
.002 w
.02381 .41938 m
.03006 .41938 L
s
P
[(4)] .01131 .41938 1 0 Mshowa
p
.002 w
.02381 .51135 m
.03006 .51135 L
s
P
[(5)] .01131 .51135 1 0 Mshowa
p
.002 w
.02381 .60332 m
.03006 .60332 L
s
P
[(6)] .01131 .60332 1 0 Mshowa
p
.001 w
.02381 .0699 m
.02756 .0699 L
s
P
p
.001 w
.02381 .08829 m
.02756 .08829 L
s
P
p
.001 w
.02381 .10668 m
.02756 .10668 L
s
P
p
.001 w
.02381 .12508 m
.02756 .12508 L
s
P
p
.001 w
.02381 .16187 m
.02756 .16187 L
s
P
p
.001 w
.02381 .18026 m
.02756 .18026 L
s
P
p
.001 w
.02381 .19865 m
.02756 .19865 L
s
P
p
.001 w
.02381 .21705 m
.02756 .21705 L
s
P
p
.001 w
.02381 .25384 m
.02756 .25384 L
s
P
p
.001 w
.02381 .27223 m
.02756 .27223 L
s
P
p
.001 w
.02381 .29062 m
.02756 .29062 L
s
P
p
.001 w
.02381 .30902 m
.02756 .30902 L
s
P
p
.001 w
.02381 .3458 m
.02756 .3458 L
s
P
p
.001 w
.02381 .3642 m
.02756 .3642 L
s
P
p
.001 w
.02381 .38259 m
.02756 .38259 L
s
P
p
.001 w
.02381 .40099 m
.02756 .40099 L
s
P
p
.001 w
.02381 .43777 m
.02756 .43777 L
s
P
p
.001 w
.02381 .45617 m
.02756 .45617 L
s
P
p
.001 w
.02381 .47456 m
.02756 .47456 L
s
P
p
.001 w
.02381 .49296 m
.02756 .49296 L
s
P
p
.001 w
.02381 .52974 m
.02756 .52974 L
s
P
p
.001 w
.02381 .54814 m
.02756 .54814 L
s
P
p
.001 w
.02381 .56653 m
.02756 .56653 L
s
P
p
.001 w
.02381 .58493 m
.02756 .58493 L
s
P
p
.001 w
.02381 .03311 m
.02756 .03311 L
s
P
p
.001 w
.02381 .01472 m
.02756 .01472 L
s
P
p
.002 w
.02381 0 m
.02381 .61803 L
s
P
P
0 0 m
1 0 L
1 .61803 L
0 .61803 L
closepath
clip
newpath
.004 w
.02381 .60332 m
.21429 .38259 L
.40476 .16187 L
.59524 .45617 L
.78571 .23544 L
.97619 .01472 L
s
% End of Graphics
MathPictureEnd
:[font = postscript; PostScript; formatAsPostScript; output; inactive; preserveAspect; pictureLeft = 34; pictureWidth = 281; pictureHeight = 173; endGroup; endGroup]
%!
%%Creator: Mathematica
%%AspectRatio: .61803
MathPictureStart
%% Graphics
/Helvetica findfont 6 scalefont setfont
% Scaling calculations
0.02381 0.095238 -0.00981 0.061313 [
[(2)] .21429 .11282 0 2 Msboxa
[(4)] .40476 .11282 0 2 Msboxa
[(6)] .59524 .11282 0 2 Msboxa
[(8)] .78571 .11282 0 2 Msboxa
[(10)] .97619 .11282 0 2 Msboxa
[(4)] .01131 .23544 1 0 Msboxa
[(6)] .01131 .35807 1 0 Msboxa
[(8)] .01131 .48069 1 0 Msboxa
[(10)] .01131 .60332 1 0 Msboxa
[ -0.001 -0.001 0 0 ]
[ 1.001 .61903 0 0 ]
] MathScale
% Start of Graphics
1 setlinecap
1 setlinejoin
newpath
[ ] 0 setdash
0 g
p
p
.002 w
.21429 .11282 m
.21429 .11907 L
s
P
[(2)] .21429 .11282 0 2 Mshowa
p
.002 w
.40476 .11282 m
.40476 .11907 L
s
P
[(4)] .40476 .11282 0 2 Mshowa
p
.002 w
.59524 .11282 m
.59524 .11907 L
s
P
[(6)] .59524 .11282 0 2 Mshowa
p
.002 w
.78571 .11282 m
.78571 .11907 L
s
P
[(8)] .78571 .11282 0 2 Mshowa
p
.002 w
.97619 .11282 m
.97619 .11907 L
s
P
[(10)] .97619 .11282 0 2 Mshowa
p
.001 w
.0619 .11282 m
.0619 .11657 L
s
P
p
.001 w
.1 .11282 m
.1 .11657 L
s
P
p
.001 w
.1381 .11282 m
.1381 .11657 L
s
P
p
.001 w
.17619 .11282 m
.17619 .11657 L
s
P
p
.001 w
.25238 .11282 m
.25238 .11657 L
s
P
p
.001 w
.29048 .11282 m
.29048 .11657 L
s
P
p
.001 w
.32857 .11282 m
.32857 .11657 L
s
P
p
.001 w
.36667 .11282 m
.36667 .11657 L
s
P
p
.001 w
.44286 .11282 m
.44286 .11657 L
s
P
p
.001 w
.48095 .11282 m
.48095 .11657 L
s
P
p
.001 w
.51905 .11282 m
.51905 .11657 L
s
P
p
.001 w
.55714 .11282 m
.55714 .11657 L
s
P
p
.001 w
.63333 .11282 m
.63333 .11657 L
s
P
p
.001 w
.67143 .11282 m
.67143 .11657 L
s
P
p
.001 w
.70952 .11282 m
.70952 .11657 L
s
P
p
.001 w
.74762 .11282 m
.74762 .11657 L
s
P
p
.001 w
.82381 .11282 m
.82381 .11657 L
s
P
p
.001 w
.8619 .11282 m
.8619 .11657 L
s
P
p
.001 w
.9 .11282 m
.9 .11657 L
s
P
p
.001 w
.9381 .11282 m
.9381 .11657 L
s
P
p
.002 w
0 .11282 m
1 .11282 L
s
P
p
.002 w
.02381 .23544 m
.03006 .23544 L
s
P
[(4)] .01131 .23544 1 0 Mshowa
p
.002 w
.02381 .35807 m
.03006 .35807 L
s
P
[(6)] .01131 .35807 1 0 Mshowa
p
.002 w
.02381 .48069 m
.03006 .48069 L
s
P
[(8)] .01131 .48069 1 0 Mshowa
p
.002 w
.02381 .60332 m
.03006 .60332 L
s
P
[(10)] .01131 .60332 1 0 Mshowa
p
.001 w
.02381 .13734 m
.02756 .13734 L
s
P
p
.001 w
.02381 .16187 m
.02756 .16187 L
s
P
p
.001 w
.02381 .18639 m
.02756 .18639 L
s
P
p
.001 w
.02381 .21092 m
.02756 .21092 L
s
P
p
.001 w
.02381 .25997 m
.02756 .25997 L
s
P
p
.001 w
.02381 .28449 m
.02756 .28449 L
s
P
p
.001 w
.02381 .30902 m
.02756 .30902 L
s
P
p
.001 w
.02381 .33354 m
.02756 .33354 L
s
P
p
.001 w
.02381 .38259 m
.02756 .38259 L
s
P
p
.001 w
.02381 .40712 m
.02756 .40712 L
s
P
p
.001 w
.02381 .43164 m
.02756 .43164 L
s
P
p
.001 w
.02381 .45617 m
.02756 .45617 L
s
P
p
.001 w
.02381 .50522 m
.02756 .50522 L
s
P
p
.001 w
.02381 .52974 m
.02756 .52974 L
s
P
p
.001 w
.02381 .55427 m
.02756 .55427 L
s
P
p
.001 w
.02381 .57879 m
.02756 .57879 L
s
P
p
.001 w
.02381 .08829 m
.02756 .08829 L
s
P
p
.001 w
.02381 .06377 m
.02756 .06377 L
s
P
p
.001 w
.02381 .03924 m
.02756 .03924 L
s
P
p
.001 w
.02381 .01472 m
.02756 .01472 L
s
P
p
.002 w
.02381 0 m
.02381 .61803 L
s
P
P
0 0 m
1 0 L
1 .61803 L
0 .61803 L
closepath
clip
newpath
.004 w
.02381 .35807 m
.11905 .21092 L
.21429 .06377 L
.30952 .25997 L
.40476 .11282 L
.5 .30902 L
.59524 .16187 L
.69048 .01472 L
.78571 .21092 L
.88095 .40712 L
.97619 .60332 L
s
% End of Graphics
MathPictureEnd
:[font = text; inactive; preserveAspect]
Lastly, we count the proportion to reach the top boundary first and the
average time to reach a boundary for this sample:
:[font = input; Cclosed; preserveAspect; startGroup]
p1 = ProportiontoHitUpper =
N[Sum[u[i], {i,1,trials}] / trials]
:[font = output; output; inactive; preserveAspect; endGroup]
0.4833333333333333
;[o]
0.483333
:[font = input; Cclosed; preserveAspect; startGroup]
SampleMeanTime =
N[Sum[t[i], {i,1,trials}] / trials]
:[font = output; output; inactive; noKeepOnOnePage; preserveAspect; endGroup]
4.766666666666667
;[o]
4.76667
:[font = text; inactive; preserveAspect]
When using h = k = 1, we can observe the accuracy of the approximations
by adding codes which compute the known real values :
:[font = input; Cclosed; preserveAspect; startGroup]
RealProb=If[p==.5,m/n,(1-(q/p)^m)/(1-(q/p)^n)]
:[font = output; output; inactive; preserveAspect; endGroup]
0.1833692373976734
;[o]
0.183369
:[font = input; Cclosed; preserveAspect; startGroup]
AveTime=If[p==.5,m(n-m),
n/(p-q)*(1-(q/p)^m)/(1-(q/p)^n) - m/(p-q)]
:[font = output; output; inactive; preserveAspect; endGroup; endGroup]
20.83153813011633
;[o]
20.8315
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Confidence Intervals
:[font = text; inactive; preserveAspect]
We can now use the data from this sample of random walks to create a
confidence interval for the probability P[m] of reaching n first when starting
at height m. The confidence interval for P[m] is centered around the sample
proportion p1. Since this measurement is a proportion, the true variance is
known to be P[m](1-P[m]) [HoTa93, p.151], which has a maximum value of
.25 since 0 ² P[m] ² 1. By taking the square root, we obtain a maximum standard
deviationof .50. The variance of the sample proportion however equals the
true variance divided by sample size [HoTa93, p. 267]; hence, the confidence
interval is centered at p1 with standard deviation .5/Sqrt[trials].
:[font = text; inactive; preserveAspect]
Using the Statistics`ConfidenceIntervals` package, we now compute a 95%
confidence interval for P[m] :
:[font = input; preserveAspect]
<c]]
:[font = output; output; inactive; dontNoPageBreakInGroup; preserveAspect; endGroup]
{0.4267540466295248, 0.5399126200371419}
;[o]
{0.426754, 0.539913}
:[font = text; inactive; preserveAspect]
Similarly, we can construct a confidence interval, centered around the
sample mean xbar, for the average time µ to reach a boundary. Though we
would now use the unbiased sample deviation sd as an estimate for the
unknown standard deviation s of the time to reach a boundary. However,
there are some other factors to consider before we do so.
:[font = text; inactive; preserveAspect]
Unlike with the sample proportion, we do not have an upper bound for
the variance of the time it takes to reach a boundary, nor do we know if this
variance is even finite. If the variance is finite, then we may apply the Central
Limit Theorem to construct a confidence interval based on a normal distribution
about the sample mean. However, we must be sure that the sample size trials
is large enough so that sd is a good estimate of s and so that the the standard
normal distribution is a good approximation of the distribution of
(xbar - µ ) / (sd/Sqrt[trials]) .
:[font = text; inactive; preserveAspect]
To test these assumptions, it is a better practice to run our simulations in
batches, say 30 batches of 100 trials each. We then consider the 30 batch means
as a new random sample. If these 30 measurements seem normally
distributed, then it is good evidence that the confidence interval constructed
will be valid. The commands below create and display such batch means.
:[font = input; preserveAspect]
Do[Do[j=0;While[0c]]
:[font = output; output; inactive; preserveAspect; endGroup; endGroup]
{4.789362890917982, 5.003303775748682}
;[o]
{4.78936, 5.0033}
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Matrix Solution
:[font = text; inactive; preserveAspect]
We now consider the special case when the upward and downward steps h
and k are integers. We shall describe an algebraic solution for the probability
P[m] of reaching height n before height 0 when starting at height m, for
0 ² m ² n. We first note that P[n] = 1, since if we start at n, then we reach n
before 0 with probability 1. Thus, we also define P[k] = 1 for k > n. Moreover,
P[0] = 0, since we cannot reach n before 0 if we start at 0, and we also define
P[k] = 0 for k < 0. So, we must solve for P[1] through P[n-1]. We shall do so by
solving n-1 equations having these n-1 unknowns. To avoid consideration of
numerous cases, we assume that k+h < n-1. Simple modifications can be made
if k+h ³ n - 1.
:[font = text; inactive; preserveAspect]
When starting at height m, for 1 ² m ² n-1, after one step we are at height
m-k with probability q or at height m+h with probability p. Then we must
reach n before 0 when starting at m-k or when starting at m+h. Hence,
:[font = text; inactive; preserveAspect]
P[m] = q P[m-k] + p P[m+h].
:[font = text; inactive; preserveAspect]
In particular, the first k equations are
:[font = text; inactive; preserveAspect]
P[1] = q P[1-k] + p P[1+h] = p P[1+h]
P[2] = q P[2-k] + p P[2+h] = p P[2+h]
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
P[k] = q P[k-k] + p P[n-1] = p P[k+h] .
:[font = text; inactive; preserveAspect]
Since k+h < n-1, the next n-1-h-k equations are
:[font = text; inactive; preserveAspect]
P[k+1] = q P[1] + p P[k+1+h]
. . . . . . . . . . . . .
. . . . . . . . . . . . .
P[n-1-h] = q P[n-1-h-k] + p P[n-1],
:[font = text; inactive; preserveAspect]
and the last h equations are
:[font = text; inactive; preserveAspect]
P[n-h] = q P[n-h-k] + p P[n] = q P[n-h-k] + p
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
P[n-1] = q P[n-1-k] + p P[n-1+h] = q P[n-1-k] + p
:[font = text; inactive; preserveAspect]
The first k equations can be written as follows:
:[font = text; inactive; preserveAspect]
P[1] - p P[1+h] = 0
P[2] - p P[2+h] = 0
. . . .
P[k] - p P[k+h] = 0
:[font = text; inactive; preserveAspect]
The next n-1-h-k equations can be written as follows:
:[font = text; inactive; preserveAspect]
-q P[1] + P[k+1] - p P[k+1+h] = 0
. . . . . . . . . .
-q P[n-1-h-k] + P[n-1-h] - p P[n-1] = 0 .
:[font = text; inactive; preserveAspect]
Finally, the last h equations can be written as
:[font = text; inactive; preserveAspect]
- q P[n-h-k] + P[n-h] = p
. . . . . .
- q P[n-1-k] + P[n-1] = p .
:[font = text; inactive; preserveAspect]
Hence, we obtain the following (n-1)x n augmented matrix which represents
the n-1 equations, having P[1] through P[n-1] as the n-1 unknowns, where the
-p is h columns to the right of the 1 in each of the first n-1-h rows, the -q is k
columns to the left of the 1 in the last n-1-k rows, and p occurs in the right-
most column in the last h rows.
:[font = input; preserveAspect]
1 0 0 0 . -p 0 . . . . 0 0 0 | 0
0 1 0 0 . 0 -p 0 . . . 0 0 0 | 0
0 0 1 0 . 0 0 -p 0 . . 0 0 0 | 0
0 0 0 1 . 0 0 0 -p 0 . 0 0 0 | 0
. . . . . . . . . . . . . . | 0
-q 0 0 . . 1 0 0 0 . -p . 0 0 | 0
0 -q 0 0 . . 1 0 0 0 . -p 0 0 | 0
0 . -q 0 0 . . 1 0 0 0 . -p 0 | 0
. . . . . . . . . . . . . . | 0
0 . . -q 0 0 . . 1 . . . 0 | p
0 . . . -q 0 0 . . 1 . . 0 | p
0 . . . . . -q 0 0 . 1 . 0 | p
0 . . . . . . -q 0 0 . . 1 . | p
0 . . . . . . -q 0 0 . . 1 | p
;[s]
1:0,0;743,-1;
1:1,10,8,Courier,0,12,0,0,0;
:[font = text; inactive; preserveAspect]
If we let A denote the (n-1)x(n-1) matrix of coefficients on the left side of
the augmented matrix, and let B denote the right-most column, then the
system can be represented by the equation AX = B, where X denotes the
column of unknowns P[1],. . ., P[n-1]. The solution for X is then given by
X = A-1 B.
;[s]
3:0,0;307,1;310,2;315,-1;
3:1,11,8,Times,0,12,0,0,0;1,11,8,Times,32,12,0,0,0;1,11,8,Times,0,12,0,0,0;
:[font = text; inactive; preserveAspect]
Of course it is not practical to solve the system of equations by hand; thus,
we can solve the system numerically for specific p and q or theoretically for
arbitrary p and q by using a short program.
:[font = text; inactive; preserveAspect]
Again, we first enter the necessary parameters and we note that this
technique will provide a solution for all starting heights m; so, we can clear
this variable. Moreover, we have changed the values of h and k to integers.
:[font = input; preserveAspect]
Clear[m,a,b];
p = .40;
q = 1-p;
StepsUp = h = 3;
StepsDown = k = 2;
UpperBoundary = n = 10;
:[font = text; inactive; preserveAspect]
The following commands create the (n-1)x(n-1) matrix of coefficients A:
:[font = input; preserveAspect]
a=IdentityMatrix[n-1];
Do[a[[i,i-k]]=-q,{i,k+1,n-1}];
Do[a[[i,i+h]]=-p,{i,1,n-1-h}];
:[font = input; Cclosed; preserveAspect; startGroup]
a
:[font = output; output; inactive; preserveAspect; endGroup]
{{1, 0, 0, -0.4, 0, 0, 0, 0, 0},
{0, 1, 0, 0, -0.4, 0, 0, 0, 0},
{-0.6, 0, 1, 0, 0, -0.4, 0, 0, 0},
{0, -0.6, 0, 1, 0, 0, -0.4, 0, 0},
{0, 0, -0.6, 0, 1, 0, 0, -0.4, 0},
{0, 0, 0, -0.6, 0, 1, 0, 0, -0.4},
{0, 0, 0, 0, -0.6, 0, 1, 0, 0},
{0, 0, 0, 0, 0, -0.6, 0, 1, 0},
{0, 0, 0, 0, 0, 0, -0.6, 0, 1}}
;[o]
{{1, 0, 0, -0.4, 0, 0, 0, 0, 0},
{0, 1, 0, 0, -0.4, 0, 0, 0, 0},
{-0.6, 0, 1, 0, 0, -0.4, 0, 0, 0},
{0, -0.6, 0, 1, 0, 0, -0.4, 0, 0},
{0, 0, -0.6, 0, 1, 0, 0, -0.4, 0},
{0, 0, 0, -0.6, 0, 1, 0, 0, -0.4},
{0, 0, 0, 0, -0.6, 0, 1, 0, 0},
{0, 0, 0, 0, 0, -0.6, 0, 1, 0},
{0, 0, 0, 0, 0, 0, -0.6, 0, 1}}
:[font = text; inactive; preserveAspect]
The next command creates the (n-1)x 1 matrix B :
:[font = input; Cclosed; preserveAspect; startGroup]
b=Table[If[i<=n-h-1,0,p],{i,1,n-1},{j,1,1}]
:[font = output; output; inactive; preserveAspect; endGroup]
{{0}, {0}, {0}, {0}, {0}, {0}, {0.4}, {0.4}, {0.4}}
;[o]
{{0}, {0}, {0}, {0}, {0}, {0}, {0.4}, {0.4}, {0.4}}
:[font = text; inactive; preserveAspect]
Lastly, we solve the system AX = B for X using the LinearSolve command.
We then display a table of values {m,P[m]}, for 1 ² m ² n-1, which gives the
probabilities of reaching n first when starting at height m. We have displayed
the output for the last command which occurs using our initial values.
;[s]
3:0,0;51,1;62,2;301,-1;
3:1,11,8,Times,0,12,0,0,0;1,10,8,Courier,1,12,0,0,0;1,11,8,Times,0,12,0,0,0;
:[font = input; preserveAspect]
Soln = LinearSolve[a, b];
:[font = input; Cclosed; preserveAspect; startGroup]
Table[{m, Soln[[m]]}, {m,1,n-1}]
:[font = output; output; inactive; preserveAspect; endGroup; endGroup]
{{1, {0.1574099955771782}}, {2, {0.1946041574524546}},
{3, {0.3193149933657673}},
{4, {0.3935249889429455}},
{5, {0.4865103936311366}},
{6, {0.5621724900486509}},
{7, {0.6919062361786819}},
{8, {0.7373034940291907}}, {9, {0.815143741707209}}}
;[o]
{{1, {0.15741}}, {2, {0.194604}}, {3, {0.319315}},
{4, {0.393525}}, {5, {0.48651}}, {6, {0.562172}},
{7, {0.691906}}, {8, {0.737303}}, {9, {0.815144}}}
:[font = section; inactive; Cclosed; preserveAspect; fontName = "Palatino"; startGroup]
Other Matrix Solutions
:[font = text; inactive; preserveAspect]
We now let R[m] be the probability of reaching (or exceeding) a boundary of 0
or n when starting at height m. With only one modification, we can use the
same matrix technique as above to verify that R[m] = 1 for all m. We note that
R[0] = R[n] = 1, since if we start at 0 or n then we reach 0 or n. Likewise R[k] = 1
for k < 0 and k > n. The only modification then is that the first k equations can
be written with the right hand side equaling q rather than 0. We need only
create a new right-most column for the augmented matrix and then solve for
R[1],. . .,R[m] :
:[font = input; Cclosed; preserveAspect; startGroup]
b2 = Table[If[i<=k,q,If[i<=n-h-1,0,p]],
{i,1,n-1}, {j,1,1}]
:[font = output; output; inactive; preserveAspect; endGroup]
{{0.6}, {0.6}, {0}, {0}, {0}, {0}, {0.4}, {0.4}, {0.4}}
;[o]
{{0.6}, {0.6}, {0}, {0}, {0}, {0}, {0.4}, {0.4}, {0.4}}
:[font = input; Cclosed; preserveAspect; startGroup]
Soln2 = LinearSolve[a, b2]
:[font = output; output; inactive; preserveAspect; endGroup]
{{1.}, {1.}, {1.}, {1.}, {1.}, {1.}, {1.}, {1.}, {1.}}
;[o]
{{1.}, {1.}, {1.}, {1.}, {1.}, {1.}, {1.}, {1.}, {1.}}
:[font = text; inactive; preserveAspect]
With another modification, we can find the average number of steps needed
to reach a boundary. If we let T[m] denote the average number when starting
at height m, then T[0] = T[n] = 0 since it does not require a step to reach the
boundary if we start at a boundary. Likewise, T[k] = 0 for k < 0 and k > n. For
0 < m < n however, we take 1 step and then we proceed from height m-k with
probability q or from height m+h with probability p. Thus,
T[m] = 1 + q T[m-k] + p T[m+h]. Since, T[k] = 0 for k ² 0 and k ³ n, we will obtain
the same system of equations as before; but now, each entry in the right-most
column will be 1. Again we need only create a new column and then solve
for each T[m].
:[font = input; Cclosed; preserveAspect; startGroup]
b3=Table[1,{i,1,n-1},{j,1,1}]
:[font = output; output; inactive; preserveAspect; endGroup]
{{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}}
;[o]
{{1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}, {1}}
:[font = input; preserveAspect]
Soln3 = LinearSolve[a, b3];
:[font = input; Cclosed; preserveAspect; startGroup]
Table[{m,Soln3[[m]]},
{m, 1, n-1}]
:[font = output; output; inactive; preserveAspect; endGroup]
{{1, {2.864449358690844}}, {2, {3.217602830605927}},
{3, {4.812674038036267}}, {4, {4.661123396727112}},
{5, {5.544007076514817}}, {6, {5.235011057054401}},
{7, {4.326404245908889}}, {8, {4.14100663423264}},
{9, {3.595842547545333}}}
;[o]
{{1, {2.86445}}, {2, {3.2176}}, {3, {4.81267}},
{4, {4.66112}}, {5, {5.54401}}, {6, {5.23501}},
{7, {4.3264}}, {8, {4.14101}}, {9, {3.59584}}}
:[font = text; inactive; preserveAspect; fontName = "Palatino"]
For h = k = 1, we could also verify the matrix solutions for P[m] and
T[m] by adding commands which list the previously known formulas:
:[font = input; preserveAspect]
Table[{m,If[p==.5,m/n,(1-(q/p)^m)/(1-(q/p)^n)]},
{m,1,n-1}]
:[font = input; preserveAspect; endGroup]
Table[{m,If[p==.5,m(n-m),
n/(p-q)*(1-(q/p)^m)/(1-(q/p)^n) - m/(p-q)]},{m,1,n-1}]
:[font = section; inactive; Cclosed; preserveAspect; fontName = "Palatino"; startGroup]
Application To Gambling
:[font = text; inactive; preserveAspect; fontName = "Palatino"]
A famous application of the ÒDrunkard's WalkÓ is the ÒGambler's RuinÓ
problem. We first find the probability Q[m] = 1 - P[m] of reaching 0 before
reaching n. This setting can now be interpreted as the probability of Ògoing
brokeÓ before reaching a goal of $n if starting with $m when one either wins
$h with probability p or loses $k with probability q on each bet. For h = k = 1,
we of course have a closed-form expression: For p = .5, Q[m] = (n-m)/n and for
p .5, Q[m] = ((q/p)^m - (q/p)^n)/(1 - (q/p)^n). Then for p ² q (or
equivalently p ² .50), lim Q[m] = 1 for all m.
n->°
:[font = text; inactive; preserveAspect]
This limit is the famous Gambler's Ruin: When one's chance of winning
on a single bet is no more than .50, then the persistent gambler will eventually
go broke with probability 1.
:[font = text; inactive; preserveAspect]
In real gambling situations though, one usually expects to win more than
$1 for each dollar bet unless the probability p of winning is close to .50. For a $1
bet with a payoff of $h, one's average earnings are -1*q + h*p. For a fair game,
the average earnings should be 0; hence, h = q/p. Since most games are in fact
unfair, so that the house can make a profit, we shall assume a payoff of $h
which equals the next integer below q/p. For p = q though, we should assume
a fair game with h = 1.
:[font = text; inactive; preserveAspect]
For example, if p = .20 and q = .80, then q/p = 4 and we take h = 3. We then
let k = 1 in the original matrix program. In other words, we either win $3 with
probability .20 or lose $1 with probability .80. Now suppose we wish to reach
$30 before going broke with probability .60. Is it possible? If so how much
money do we initially need? By running the program with p = .20, h = 3, k = 1,
and n = 30, we obtain the following final output :
:[font = output; output; inactive; preserveAspect]
{{1, {0.002054711548168111571}},
{2, {0.004419499844303345578}},
{3, {0.007141155905222170139}},
{4, {0.01027355774084055785}},
{5, {0.01387865302884428161}},
{6, {0.01802778014889746838}},
{7, {0.0228031650833141087}},
{8, {0.02829903418085917664}},
{9, {0.03462428862911021547}},
{10, {0.04190470482098066995}},
{11, {0.05028251057103944841}},
{12, {0.05992530642211437079}},
{13, {0.0710263695884624879}},
{14, {0.08379373357127456221}},
{15, {0.09849648982641406034}},
{16, {0.1154306222538549563}},
{17, {0.1348631895025228594}},
{18, {0.1573075148469720528}},
{19, {0.1831671519636185402}},
{20, {0.2125934584971944719}},
{21, {0.2470848162247688264}},
{22, {0.2866057004302044899}},
{23, {0.3302986846314981986}},
{24, {0.3850502471350662443}},
{25, {0.4446892372519471438}},
{26, {0.5050706214366730337}},
{27, {0.604056497149338427}},
{28, {0.6832451977194707417}},
{29, {0.7465961581755765933}}}
;[o]
{{1, {0.00205471}}, {2, {0.0044195}}, {3, {0.00714116}},
{4, {0.0102736}}, {5, {0.0138787}}, {6, {0.0180278}},
{7, {0.0228032}}, {8, {0.028299}}, {9, {0.0346243}},
{10, {0.0419047}}, {11, {0.0502825}}, {12, {0.0599253}},
{13, {0.0710264}}, {14, {0.0837937}}, {15, {0.0984965}},
{16, {0.115431}}, {17, {0.134863}}, {18, {0.157308}},
{19, {0.183167}}, {20, {0.212593}}, {21, {0.247085}},
{22, {0.286606}}, {23, {0.330299}}, {24, {0.38505}},
{25, {0.444689}}, {26, {0.505071}}, {27, {0.604056}},
{28, {0.683245}}, {29, {0.746596}}}
:[font = text; inactive; preserveAspect; fontName = "Palatino"; endGroup]
Hence, if we start with $27, we have just slightly more than a 60% chance
of reaching $30 before going broke. By experimenting with larger and larger
values of n, it becomes evident that there is an upper bound to the probability
of reaching $n before going broke. With the above figures, this upper bound
appears to be .75. Moreover, we cannot contradict the Gambler's Ruin. That is,
for a fixed initial amount m, if n-> °, then P[m]-> 0 and hence Q[m] -> 1. One
will still eventually go broke if one persists in gambling. However, if one
decides to Òquit while you're aheadÓ, one can determine the precise probability
of reaching the goal of $n before going broke!
:[font = section; inactive; Cclosed; preserveAspect; fontName = "Palatino"; startGroup]
Exercise
:[font = text; inactive; preserveAspect; fontName = "Palatino"; endGroup]
Find an estimation and confidence interval for the average Òarea under the
curveÓ A[m] when starting at height m and stopping upon reaching height n
or height 0. Adjust the matrix solution to find the exact value of A[m] when
the steps h and k are integers.
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
Acknowledgments
:[font = text; inactive; preserveAspect; endGroup]
The author would like to thank the referees for their suggestions and
pointers which made the programming more efficient and elegant.
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
References
:[font = text; inactive; preserveAspect]
[Chu75] K.L. Chung. Elementary Probability Theory with Stochastic
Processes . Springer-Verlag, New York, 1975.
:[font = text; inactive; preserveAspect]
[Fel57] W. Feller. An Introduction to Probability Theory and Its Applications
Volume I, Second Edition. John Wiley & Sons, Tokyo, 1957.
:[font = text; inactive; preserveAspect; endGroup]
[HoTa93] R.V. Hogg and E.A. Tanis. Probability and Statistical Inference
Fourth Edition. Macmillan Publishing Co., New York, 1993.
:[font = section; inactive; Cclosed; preserveAspect; startGroup]
About the Author
:[font = text; inactive; preserveAspect]
David K. Neal is an associate professor of mathematics at Western Kentucky
University where he has taught for six years. His mathematical interests are
in stochastic analysis and general probability. In recent years, he has enjoyed
enriching his courses with a variety computer projects using Mathematica .
His other hobbies include music, sports, and felines.
:[font = text; inactive; preserveAspect; endGroup; endGroup]
David K. Neal
Department of Mathematics
Western Kentucky University
Bowling Green, KY 42101
nealdk @wkuvx1.wku.edu
^*)