ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sat, 11 Nov 2017 10:47:52 +0100Sage Probabilityhttps://ask.sagemath.org/question/10584/sage-probability/Please help me write a sage function for the following problem:
Suppose a gambler starts with R1000. He repeatedly plays a game where the initial bet is R100. The gambler has a 60% chance of losing the game (in which case he gets nothing), but a 40% chance of winning and doubling his money (ie, he gets back R200). If the gambler loses all of his money, he is ruined and stops playing. If the gambler reaches R2000, he also stops playing. Otherwise, the gambler keeps playing. What is the probability that the gambler stops with R2000? What is the probability he is ruined?
Write a simulation to empirically determine the probabilities. Run your simulation 10000 times. What are the probabilities for ruin and winning?Tue, 01 Oct 2013 12:42:12 +0200https://ask.sagemath.org/question/10584/sage-probability/Answer by dan_fulea for <p>Please help me write a sage function for the following problem:</p>
<p>Suppose a gambler starts with R1000. He repeatedly plays a game where the initial bet is R100. The gambler has a 60% chance of losing the game (in which case he gets nothing), but a 40% chance of winning and doubling his money (ie, he gets back R200). If the gambler loses all of his money, he is ruined and stops playing. If the gambler reaches R2000, he also stops playing. Otherwise, the gambler keeps playing. What is the probability that the gambler stops with R2000? What is the probability he is ruined?</p>
<p>Write a simulation to empirically determine the probabilities. Run your simulation 10000 times. What are the probabilities for ruin and winning?</p>
https://ask.sagemath.org/question/10584/sage-probability/?answer=38085#post-id-38085
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.
Mon, 26 Jun 2017 03:44:24 +0200https://ask.sagemath.org/question/10584/sage-probability/?answer=38085#post-id-38085Comment by tmonteil for <p>OK, let us do the homework.</p>
<p>The simulation first, as it was required.</p>
<pre><code>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 )
</code></pre>
<p>Results after running the code three times:</p>
<pre><code>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
</code></pre>
<p>OK, this was the homework.</p>
<p>Theoretical check.
The probability model is a tree, i cannot paint it here,
but in a plain ascii world it looks like...</p>
<pre><code> * +5
/
*
/ \
* * +3
/ \ /
* *
/ \ / \
* * * +1
/ \ / \ /
* * *
\ / \ / \
* * * -1
\ / \ /
* *
\ / \
* * -3
\ /
*
\
* -5
</code></pre>
<p>The up-moves have probability $p=\frac 25$.</p>
<p>The down-moves have probability $q=1-p=\frac 35$.</p>
<p>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$.</p>
<p>We need these two probabilities.
(And implicitly the probability that the gambler lives his live in this
one game.)</p>
<p>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</p>
<pre><code> *--*--*--...--*--*--*--...--*--*--*
-10 -9 -8 -1 0 +1 +8 +9 +10
</code></pre>
<p>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$.</p>
<p>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.)</p>
<p>The problem can now be solved theoretically. The associated Markov matrix is</p>
<pre><code>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()
</code></pre>
<p>This gives:</p>
<pre><code>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]
</code></pre>
<p>(To force the explicit print, one has to know and use the <code>str</code> method.)</p>
<p>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).</p>
<p>So we can simply ask sage for...</p>
<pre><code>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]
</code></pre>
<p>This gives:</p>
<pre><code>Probability of win is numerically:
w[20] ~ 0.017045927454929835366970186273367402993025152730844139630117
Probability of ruin is numerically:
w[ 0] ~ 0.98295407254507016463302981372663259700697484726915586036990
</code></pre>
<p>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:</p>
<pre><code>sage: c = 3/2
sage: x0 = (c^10-1)/(c^20-1)
sage: x0
1024/60073
sage: RealField( 200 )( x0 )
0.017045927454929835366970186273367402993025152730844139630117
</code></pre>
<p>and coincides with the numerically computed value.</p>
<p>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.</p>
https://ask.sagemath.org/question/10584/sage-probability/?comment=39496#post-id-39496It is a very good idea to solve homeworks 3+ years after they have been asked. BY, the, way, it is also interesting to compare the speed of convergence between repeated simulations and matrix exponentiation.Sat, 11 Nov 2017 10:47:50 +0100https://ask.sagemath.org/question/10584/sage-probability/?comment=39496#post-id-39496Comment by tmonteil for <p>OK, let us do the homework.</p>
<p>The simulation first, as it was required.</p>
<pre><code>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 )
</code></pre>
<p>Results after running the code three times:</p>
<pre><code>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
</code></pre>
<p>OK, this was the homework.</p>
<p>Theoretical check.
The probability model is a tree, i cannot paint it here,
but in a plain ascii world it looks like...</p>
<pre><code> * +5
/
*
/ \
* * +3
/ \ /
* *
/ \ / \
* * * +1
/ \ / \ /
* * *
\ / \ / \
* * * -1
\ / \ /
* *
\ / \
* * -3
\ /
*
\
* -5
</code></pre>
<p>The up-moves have probability $p=\frac 25$.</p>
<p>The down-moves have probability $q=1-p=\frac 35$.</p>
<p>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$.</p>
<p>We need these two probabilities.
(And implicitly the probability that the gambler lives his live in this
one game.)</p>
<p>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</p>
<pre><code> *--*--*--...--*--*--*--...--*--*--*
-10 -9 -8 -1 0 +1 +8 +9 +10
</code></pre>
<p>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$.</p>
<p>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.)</p>
<p>The problem can now be solved theoretically. The associated Markov matrix is</p>
<pre><code>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()
</code></pre>
<p>This gives:</p>
<pre><code>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]
</code></pre>
<p>(To force the explicit print, one has to know and use the <code>str</code> method.)</p>
<p>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).</p>
<p>So we can simply ask sage for...</p>
<pre><code>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]
</code></pre>
<p>This gives:</p>
<pre><code>Probability of win is numerically:
w[20] ~ 0.017045927454929835366970186273367402993025152730844139630117
Probability of ruin is numerically:
w[ 0] ~ 0.98295407254507016463302981372663259700697484726915586036990
</code></pre>
<p>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:</p>
<pre><code>sage: c = 3/2
sage: x0 = (c^10-1)/(c^20-1)
sage: x0
1024/60073
sage: RealField( 200 )( x0 )
0.017045927454929835366970186273367402993025152730844139630117
</code></pre>
<p>and coincides with the numerically computed value.</p>
<p>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.</p>
https://ask.sagemath.org/question/10584/sage-probability/?comment=39497#post-id-39497It is a very good idea to solve homeworks 3+ years after they have been asked. BYy the, way, it is also interesting to compare the speed of convergence between repeated simulations and matrix exponentiation.Sat, 11 Nov 2017 10:47:52 +0100https://ask.sagemath.org/question/10584/sage-probability/?comment=39497#post-id-39497