Ask Your Question
0

Sage Probability

asked 2013-10-01 12:42:12 +0100

Komolo gravatar image

updated 2015-05-24 15:28:00 +0100

FrédéricC gravatar image

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?

edit retag flag offensive close merge delete

1 Answer

Sort by » oldest newest most voted
4

answered 2017-06-26 03:44:24 +0100

dan_fulea gravatar image

updated 2017-06-26 03:45:47 +0100

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.

edit flag offensive delete link more

Comments

It 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.

tmonteil gravatar imagetmonteil ( 2017-11-11 10:47:50 +0100 )edit

It 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.

tmonteil gravatar imagetmonteil ( 2017-11-11 10:47:52 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2013-10-01 12:42:12 +0100

Seen: 758 times

Last updated: Jun 26 '17