OK, let us do the homework.
The simulation first, as it was required.
import random
def trial( p=0.4 ):
"""
repeat a binomially distributed game with parameters
p = 2/5 (default), q = 3/5
till the player either wins a total of 10
or loses a total of ten.
At each step,
winning has probability p and the win +1,
losing has probability q and the "win" -1,
q = 1-p.
"""
win = 0
while abs(win) < 10:
if random.uniform( 0, 1 ) < p:
win += 1
else:
win -= 1
if win == 10:
return 1 # win
return 0 # ruin
# let us check the win rate...
# we implicitly assume that the game will be surely over
wins = 0
N = 10000
for tryCounter in xrange( N ):
wins += trial()
print "Number or wins this time = %s" % wins
print "Winning rate this time = %s" % float( wins / N )
Results after running the code three times:
Number or wins this time = 163
Winning rate this time = 0.0163
Number or wins this time = 184
Winning rate this time = 0.0184
Number or wins this time = 159
Winning rate this time = 0.0159
OK, this was the homework.
Theoretical check.
The probability model is a tree, i cannot paint it here,
but in a plain ascii world it looks like...
* +5
/
*
/ \
* * +3
/ \ /
* *
/ \ / \
* * * +1
/ \ / \ /
* * *
\ / \ / \
* * * -1
\ / \ /
* *
\ / \
* * -3
\ /
*
\
* -5
The up-moves have probability $p=\frac 25$.
The down-moves have probability $q=1-p=\frac 35$.
At each step there is a win of 1 or -1. (Instead of 100.)
Only a portion of the combinatorial model is shown.
The gambler wins and stops when it reaches $+10$.
The gambler is loosing and the game stops when it reaches $-10$.
We need these two probabilities.
(And implicitly the probability that the gambler lives his live in this
one game.)
So the binomial model has to be truncated "horizontally" when it hits $\pm 10$.
In fact we change this image, better use the Markov model based on
*--*--*--...--*--*--*--...--*--*--*
-10 -9 -8 -1 0 +1 +8 +9 +10
and the states $\pm 10$ are "terminal",
and from each other state there is only a one unit step move
possible, to the right with probability $p=\frac 25$,
to the left with the probability $q=1-p=\frac 35$.
It is natural to make the states $\pm 10$ absorbant, and put
the probability one to stay in there, after reaching the one or the other
final state. (Convention, in the sense of the statement.)
The problem can now be solved theoretically. The associated Markov matrix is
var( 'p,q' );
def MarkovMatrix( p,q, rows=21, ring=SR ):
A = matrix( ring, rows, rows )
A[0,0] = 1
A[rows-1,rows-1] = 1
for k in range( 1, rows-1 ):
A[k,k-1] = q
A[k,k+1] = p
return A
A = MarkovMatrix( p,q )
print "The matrix A is:"
print A
print A.str()
This gives:
The matrix A is:
21 x 21 dense matrix over Symbolic Ring
[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[q 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 q 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 q 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 q 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 q 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 q 0 p 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 q 0 p 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 q 0 p 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 q 0 p 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 q 0 p 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 q 0 p 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 q 0 p 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 q 0 p 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 q 0 p 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 q 0 p 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 q 0 p 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 q 0 p 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 q 0 p 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 q 0 p]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]
(To force the explicit print, one has to know and use the str method.)
The exercise wants from us the value of
$$ w = \lim_{n\to\infty} v A^n = v \lim_{n\to\infty} A^n \ ,$$
where $v$ is the initial distribution of probability,
the row vector with entries $0,0,\dots,0,1,0,\dots,0,0$
where the one probability weight corresponds to the state $0$
(in the middle of the sketched chain).
So we can simply ask sage for...
R = RealField( 200 )
v = vector( R, 21 ); v[10] = 1
A = MarkovMatrix( 0.4, 0.6, rows=21, ring=R )
w = v*A^(2^20)
print "Probability of win is numerically:\nw[20] ~ %s" % w[20]
print "Probability of ruin is numerically:\nw[ 0] ~ %s" % w[ 0]
This gives:
Probability of win is numerically:
w[20] ~ 0.017045927454929835366970186273367402993025152730844139630117
Probability of ruin is numerically:
w[ 0] ~ 0.98295407254507016463302981372663259700697484726915586036990
We finally compute the probability of win exactly.
Let $x(k)$, $-10\le k\le +10$ be the probability to win, starting from the
state $k$. Then we have of course $x(-10)=0$, $x(+10)=1$, and for the
other possible values of the state $k$ we have the one step transition
equation:
$$ x(k) = px(k+1)+qx(k-1)\ .$$
Solving this system of equations we get the solution.
Maybe it is better to prefer to consider the recursive sequence defined by the above
recursion. The characteristic equation associated is $X = pX^2+q$, it has the
roots $1$ and $c:=q/p=\frac 32$ in our case. The general term is then given
by a formula of the shape
$$ x(n)= A\cdot 1 ^n + Bc^n=A+Bc^n\ .$$
We determine $A,B$, since we know the
values for $k=\pm10$. From
$$ A+Bc^{-10}=0\ ,$$
$$ A+Bc^{+10}=1\ ,$$
we get $B$ after subtracting, then $A$,
$$ B=\frac 1{c^{10}-c^{-10}}\ ,\ A = -c^{-10}B\ ,$$
so the needed probability is
$$ x(0)= A+Bc^0 = A+B = \frac {1-c^{-10}}{c^{10}-c^{-10}}
=\frac {c^{10}-1}{c^{20}-1}\ .$$
This is:
sage: c = 3/2
sage: x0 = (c^10-1)/(c^20-1)
sage: x0
1024/60073
sage: RealField( 200 )( x0 )
0.017045927454929835366970186273367402993025152730844139630117
and coincides with the numerically computed value.
Note: It is didactically important to provide simple code
and simple triple checks as above for
such problems, this shows how easy sage can be used and how it works in
harmony with mathematics.
