In [3]:
# generate an 11 by 11 transition matrix for a simple dice game: roll and sum until a score of 10
# We have states {0,1,2,...,9} corresponding to our scores and a state 10 corresponding to a score of 10 or more.
A=zero_matrix(QQ,11);
for i in range(6):
A[0,i+1]=1/6
A[10,10]=1;
for i in range(1,10):
for j in range(6):
roll=j+1
put = min(10,i+roll)
A[i,put]+= 1/6
A
Out[3]:
In [9]:
# output the matrix in LaTeX format
latex(A)
Out[9]:
In [12]:
# factor out a 1/6 for simpler presentation
latex((6*A))
Out[12]:
In [21]:
# find probability that the game has reached the absorbing state in 5 or fewer rolls
(A^5)[0][10]
Out[21]:
In [5]:
# work out the expected number of rolls until absorption
# A is already in canonical form, so we get the matrix Q by removing the last row and last column
Q = A[:10,:10]
Q
Out[5]:
In [13]:
N=(matrix.identity(10)-Q).inverse()
In [15]:
N
Out[15]:
In [16]:
sum(N[0][i] for i in range(10)) # note we get an exact value here
Out[16]:
In [17]:
sum(N[0][i] for i in range(10))*1. # it's always nice to report a decimal equivalent in addition to the exact value
Out[17]:
In [26]:
# lets investigate the expected number of rolls as the game limit changesA=zero_matrix(QQ,11);
R=[]
for n in range(7,50):
A=zero_matrix(QQ,n+1);
for i in range(6):
A[0,i+1]=1/6
A[n,n]=1;
for i in range(1,n):
for j in range(6):
roll=j+1
put = min(n,i+roll)
A[i,put]+= 1/6
Q = A[:n,:n]
N=(matrix.identity(n)-Q).inverse()
s=sum(N[0][i]for i in range(n))
print(n," ",s*1.)
R.append(s*1.)
list_plot(R) ## this plot is quite linear (we can prove this)! Your plot may not be so linear.
Out[26]:
In [ ]: