### load some libraries
library(kableExtra)
library(expm)      # needed for matrix power operator %^%

Chapter 11.2, Exercise 13(b), Page 424 (part of assigned homework)

Smith is in jail and has 1 dollar; he can get out on bail if he has 8 dollars.
A guard agrees to make a series of bets with him.

If Smith bets A dollars,

he wins A dollars with probability .4 ,
and loses A dollars with probability .6 .

Find the probability that he wins 8 dollars before losing all of his money if:

(b) bold strategy

he bets, each time, as much as possible but not more than necessary to bring his fortune up to 8 dollars

Here are three different ways to obtain the answer:

(b1) only three plays

Under the Bold strategy, the prisoner bets his entire bankroll each time (as it’s not possible to reach any value in \(\{5,6,7\}\) when starting from 1 dollar.)

Thus, to win the game, he must win three times:

Bet 1: 1 dollar increases to 2 dollars
Bet 2: 2 dollars increases to 4 dollars
Bet 3: 4 dollars increases to 8 dollars

A loss on any of the above bets will “wipe out” his bankroll, forcing him to stop playing.

Because his chance of winning on each bet is 0.4, his chance of winning three times is \((0.4)^3 = 0.064\) .

Thus, his probability of winning under this strategy is 6.4 percent.

(b2) Absorbing Markov Chains:

# Probabilities of win or lose
p=0.4
q=1-p

# bold Markov transition probabilities
nodes=0:8
boldmarkov=matrix(
  c(
1,0,0,0,0,0,0,0,0,   # have zero - game lost
q,0,p,0,0,0,0,0,0,   # have one -- bet one -- move to zero or two
q,0,0,0,p,0,0,0,0,   # have two -- bet two -- move to zero or four
q,0,0,0,0,0,p,0,0,   # have three -- bet THREE -- move to ZERO or SIX
q,0,0,0,0,0,0,0,p,   # have four -- bet four -- move to zero or eight
0,0,q,0,0,0,0,0,p,   # have five -- bet three -- move to two or eight
0,0,0,0,q,0,0,0,p,   # have six -- bet two -- move to four or eight
0,0,0,0,0,0,q,0,p,   # have seven -- bet one -- move to six or eight
0,0,0,0,0,0,0,0,1),  # have eight -- game won
nrow = 9, ncol = 9, byrow = T, dimnames = list(nodes,nodes)
)
# Canonical order, P
canonicalorder=c("1","2","3","4","5","6","7","0","8")
Pbold=boldmarkov[canonicalorder,canonicalorder]
Pbold

##   1   2 3   4 5   6 7   0   8
## 1 0 0.4 0 0.0 0 0.0 0 0.6 0.0
## 2 0 0.0 0 0.4 0 0.0 0 0.6 0.0
## 3 0 0.0 0 0.0 0 0.4 0 0.6 0.0
## 4 0 0.0 0 0.0 0 0.0 0 0.6 0.4
## 5 0 0.6 0 0.0 0 0.0 0 0.0 0.4
## 6 0 0.0 0 0.6 0 0.0 0 0.0 0.4
## 7 0 0.0 0 0.0 0 0.6 0 0.0 0.4
## 0 0 0.0 0 0.0 0 0.0 0 1.0 0.0
## 8 0 0.0 0 0.0 0 0.0 0 0.0 1.0

# Canonical transitions, Q
Qbold = Pbold[1:7,1:7]
Qbold

##   1   2 3   4 5   6 7
## 1 0 0.4 0 0.0 0 0.0 0
## 2 0 0.0 0 0.4 0 0.0 0
## 3 0 0.0 0 0.0 0 0.4 0
## 4 0 0.0 0 0.0 0 0.0 0
## 5 0 0.6 0 0.0 0 0.0 0
## 6 0 0.0 0 0.6 0 0.0 0
## 7 0 0.0 0 0.0 0 0.6 0

# Transient-to-absorbing, R
Rbold = Pbold[1:7,8:9]
Rbold

##     0   8
## 1 0.6 0.0
## 2 0.6 0.0
## 3 0.6 0.0
## 4 0.6 0.4
## 5 0.0 0.4
## 6 0.0 0.4
## 7 0.0 0.4

# Identity matrix, I
I = diag(7)
I

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    1    0    0    0    0    0    0
## [2,]    0    1    0    0    0    0    0
## [3,]    0    0    1    0    0    0    0
## [4,]    0    0    0    1    0    0    0
## [5,]    0    0    0    0    1    0    0
## [6,]    0    0    0    0    0    1    0
## [7,]    0    0    0    0    0    0    1

# Identity minus Q
ImQbold = I-Qbold
ImQbold

##   1    2 3    4 5    6 7
## 1 1 -0.4 0  0.0 0  0.0 0
## 2 0  1.0 0 -0.4 0  0.0 0
## 3 0  0.0 1  0.0 0 -0.4 0
## 4 0  0.0 0  1.0 0  0.0 0
## 5 0 -0.6 0  0.0 1  0.0 0
## 6 0  0.0 0 -0.6 0  1.0 0
## 7 0  0.0 0  0.0 0 -0.6 1

# N, the Fundamental Matrix
Nbold=solve(ImQbold)
Nbold

##   1   2 3    4 5   6 7
## 1 1 0.4 0 0.16 0 0.0 0
## 2 0 1.0 0 0.40 0 0.0 0
## 3 0 0.0 1 0.24 0 0.4 0
## 4 0 0.0 0 1.00 0 0.0 0
## 5 0 0.6 0 0.24 1 0.0 0
## 6 0 0.0 0 0.60 0 1.0 0
## 7 0 0.0 0 0.36 0 0.6 1

# Expected time to absorption
c = rep(1,7)
Ncbold = Nbold %*% c
Ncbold

##   [,1]
## 1 1.56
## 2 1.40
## 3 1.64
## 4 1.00
## 5 1.84
## 6 1.60
## 7 1.96

# Absorption probabilities, B
Bbold = Nbold %*% Rbold
Bbold

##       0     8
## 1 0.936 0.064
## 2 0.840 0.160
## 3 0.744 0.256
## 4 0.600 0.400
## 5 0.504 0.496
## 6 0.360 0.640
## 7 0.216 0.784

# Probability of success, starting from 1 dollar:
p1bold = Bbold["1","8"]
p1bold

## [1] 0.064

Again, the probability is 0.064 .

(b3) Repeated matrix multiplications:

p=0.4
q=1-p
nodes=0:8
boldmarkov

##     0 1   2 3   4 5   6 7   8
## 0 1.0 0 0.0 0 0.0 0 0.0 0 0.0
## 1 0.6 0 0.4 0 0.0 0 0.0 0 0.0
## 2 0.6 0 0.0 0 0.4 0 0.0 0 0.0
## 3 0.6 0 0.0 0 0.0 0 0.4 0 0.0
## 4 0.6 0 0.0 0 0.0 0 0.0 0 0.4
## 5 0.0 0 0.6 0 0.0 0 0.0 0 0.4
## 6 0.0 0 0.0 0 0.6 0 0.0 0 0.4
## 7 0.0 0 0.0 0 0.0 0 0.6 0 0.4
## 8 0.0 0 0.0 0 0.0 0 0.0 0 1.0

boldmarkov

##     0 1   2 3   4 5   6 7   8
## 0 1.0 0 0.0 0 0.0 0 0.0 0 0.0
## 1 0.6 0 0.4 0 0.0 0 0.0 0 0.0
## 2 0.6 0 0.0 0 0.4 0 0.0 0 0.0
## 3 0.6 0 0.0 0 0.0 0 0.4 0 0.0
## 4 0.6 0 0.0 0 0.0 0 0.0 0 0.4
## 5 0.0 0 0.6 0 0.0 0 0.0 0 0.4
## 6 0.0 0 0.0 0 0.6 0 0.0 0 0.4
## 7 0.0 0 0.0 0 0.0 0 0.6 0 0.4
## 8 0.0 0 0.0 0 0.0 0 0.0 0 1.0

boldresult = vector()
boldtransient = vector()
epsilon = 1e-8
for (i in 1:10) {
  boldpower = boldmarkov %^% i

  # boldmarkov raised to power i
  print(paste("i: ",i))
  print(round(boldpower,8))
  
  boldresult[i]=boldpower["1","8"]
  boldtransient[i] = sum(colSums(boldpower)[2:7])
  if (i>1)
    bolddiff=boldresult[i]-boldresult[i-1] 
  else bolddiff=999

  print(paste(i,round(boldtransient[i],10),round(boldresult[i],   8),round(bolddiff,10)))

  if (boldtransient[i] < epsilon) {
    print(paste("Absorbed at power ", i))
    break
  }

}

## [1] "i:  1"
##     0 1   2 3   4 5   6 7   8
## 0 1.0 0 0.0 0 0.0 0 0.0 0 0.0
## 1 0.6 0 0.4 0 0.0 0 0.0 0 0.0
## 2 0.6 0 0.0 0 0.4 0 0.0 0 0.0
## 3 0.6 0 0.0 0 0.0 0 0.4 0 0.0
## 4 0.6 0 0.0 0 0.0 0 0.0 0 0.4
## 5 0.0 0 0.6 0 0.0 0 0.0 0 0.4
## 6 0.0 0 0.0 0 0.6 0 0.0 0 0.4
## 7 0.0 0 0.0 0 0.0 0 0.6 0 0.4
## 8 0.0 0 0.0 0 0.0 0 0.0 0 1.0
## [1] "1 3 0 999"
## [1] "i:  2"
##      0 1 2 3    4 5 6 7    8
## 0 1.00 0 0 0 0.00 0 0 0 0.00
## 1 0.84 0 0 0 0.16 0 0 0 0.00
## 2 0.84 0 0 0 0.00 0 0 0 0.16
## 3 0.60 0 0 0 0.24 0 0 0 0.16
## 4 0.60 0 0 0 0.00 0 0 0 0.40
## 5 0.36 0 0 0 0.24 0 0 0 0.40
## 6 0.36 0 0 0 0.00 0 0 0 0.64
## 7 0.00 0 0 0 0.36 0 0 0 0.64
## 8 0.00 0 0 0 0.00 0 0 0 1.00
## [1] "2 1 0 0"
## [1] "i:  3"
##       0 1 2 3 4 5 6 7     8
## 0 1.000 0 0 0 0 0 0 0 0.000
## 1 0.936 0 0 0 0 0 0 0 0.064
## 2 0.840 0 0 0 0 0 0 0 0.160
## 3 0.744 0 0 0 0 0 0 0 0.256
## 4 0.600 0 0 0 0 0 0 0 0.400
## 5 0.504 0 0 0 0 0 0 0 0.496
## 6 0.360 0 0 0 0 0 0 0 0.640
## 7 0.216 0 0 0 0 0 0 0 0.784
## 8 0.000 0 0 0 0 0 0 0 1.000
## [1] "3 0 0.064 0.064"
## [1] "Absorbed at power  3"

# Converged power matrix
boldpower

##       0 1 2 3 4 5 6 7     8
## 0 1.000 0 0 0 0 0 0 0 0.000
## 1 0.936 0 0 0 0 0 0 0 0.064
## 2 0.840 0 0 0 0 0 0 0 0.160
## 3 0.744 0 0 0 0 0 0 0 0.256
## 4 0.600 0 0 0 0 0 0 0 0.400
## 5 0.504 0 0 0 0 0 0 0 0.496
## 6 0.360 0 0 0 0 0 0 0 0.640
## 7 0.216 0 0 0 0 0 0 0 0.784
## 8 0.000 0 0 0 0 0 0 0 1.000

# Sums of columns
print("Column sums: ")

## [1] "Column sums: "

print(round(colSums(boldpower),8))

##   0   1   2   3   4   5   6   7   8 
## 5.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.8

# bold result
print(paste("bold result: ",round(boldresult[i],8)))

## [1] "bold result:  0.064"

The probability of winning with the bold strategy is 0.064 .

Chapter 11.2, Exercise 14, Page 424

With the situation in Exercise 13, consider the strategy such that

for \(i < 4\), Smith bets \(min(i, 4-i)\), and
for \(i \ge 4\), he bets according to the bold strategy, where \(i\) is his current fortune.

Find the probability that he gets out of jail using this strategy.

(a) Absorbing Markov Chains (section 11.2):

# Probabilities of win or lose
p=0.4
q=1-p

# Markov transition probabilities
mixedmarkov=matrix(c(
1,0,0,0,0,0,0,0,0,   # have zero - game lost
q,0,p,0,0,0,0,0,0,   # have one -- bet one -- move to zero or two
q,0,0,0,p,0,0,0,0,   # have two -- bet two -- move to zero or four
0,0,q,0,p,0,0,0,0,   # have three -- bet one -- move to two or four
q,0,0,0,0,0,0,0,p,   # have four -- bet four -- move to zero or eight
0,0,q,0,0,0,0,0,p,   # have five -- bet three -- move to two or eight
0,0,0,0,q,0,0,0,p,   # have six -- bet two -- move to four or eight
0,0,0,0,0,0,q,0,p,   # have seven -- bet one -- move to six or eight
0,0,0,0,0,0,0,0,1),  # have eight -- game won
nrow = 9, ncol = 9, 
byrow = T, dimnames = list(0:8,0:8)
)
mixedmarkov

##     0 1   2 3   4 5   6 7   8
## 0 1.0 0 0.0 0 0.0 0 0.0 0 0.0
## 1 0.6 0 0.4 0 0.0 0 0.0 0 0.0
## 2 0.6 0 0.0 0 0.4 0 0.0 0 0.0
## 3 0.0 0 0.6 0 0.4 0 0.0 0 0.0
## 4 0.6 0 0.0 0 0.0 0 0.0 0 0.4
## 5 0.0 0 0.6 0 0.0 0 0.0 0 0.4
## 6 0.0 0 0.0 0 0.6 0 0.0 0 0.4
## 7 0.0 0 0.0 0 0.0 0 0.6 0 0.4
## 8 0.0 0 0.0 0 0.0 0 0.0 0 1.0

# Canonical order, P
canonicalorder=c("1","2","3","4","5","6","7","0","8")
Pmixed=mixedmarkov[canonicalorder,canonicalorder]
Pmixed

##   1   2 3   4 5   6 7   0   8
## 1 0 0.4 0 0.0 0 0.0 0 0.6 0.0
## 2 0 0.0 0 0.4 0 0.0 0 0.6 0.0
## 3 0 0.6 0 0.4 0 0.0 0 0.0 0.0
## 4 0 0.0 0 0.0 0 0.0 0 0.6 0.4
## 5 0 0.6 0 0.0 0 0.0 0 0.0 0.4
## 6 0 0.0 0 0.6 0 0.0 0 0.0 0.4
## 7 0 0.0 0 0.0 0 0.6 0 0.0 0.4
## 0 0 0.0 0 0.0 0 0.0 0 1.0 0.0
## 8 0 0.0 0 0.0 0 0.0 0 0.0 1.0

# Canonical transitions, Q
Qmixed = Pmixed[1:7,1:7]
Qmixed

##   1   2 3   4 5   6 7
## 1 0 0.4 0 0.0 0 0.0 0
## 2 0 0.0 0 0.4 0 0.0 0
## 3 0 0.6 0 0.4 0 0.0 0
## 4 0 0.0 0 0.0 0 0.0 0
## 5 0 0.6 0 0.0 0 0.0 0
## 6 0 0.0 0 0.6 0 0.0 0
## 7 0 0.0 0 0.0 0 0.6 0

# Transient-to-absorbing, R
Rmixed = Pmixed[1:7,8:9]
Rmixed

##     0   8
## 1 0.6 0.0
## 2 0.6 0.0
## 3 0.0 0.0
## 4 0.6 0.4
## 5 0.0 0.4
## 6 0.0 0.4
## 7 0.0 0.4

# Identity matrix, I
I = diag(7)
I

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    1    0    0    0    0    0    0
## [2,]    0    1    0    0    0    0    0
## [3,]    0    0    1    0    0    0    0
## [4,]    0    0    0    1    0    0    0
## [5,]    0    0    0    0    1    0    0
## [6,]    0    0    0    0    0    1    0
## [7,]    0    0    0    0    0    0    1

# Identity minus Q
ImQmixed = I-Qmixed
ImQmixed

##   1    2 3    4 5    6 7
## 1 1 -0.4 0  0.0 0  0.0 0
## 2 0  1.0 0 -0.4 0  0.0 0
## 3 0 -0.6 1 -0.4 0  0.0 0
## 4 0  0.0 0  1.0 0  0.0 0
## 5 0 -0.6 0  0.0 1  0.0 0
## 6 0  0.0 0 -0.6 0  1.0 0
## 7 0  0.0 0  0.0 0 -0.6 1

# N, the Fundamental Matrix
Nmixed=solve(ImQmixed)
Nmixed

##   1   2 3    4 5   6 7
## 1 1 0.4 0 0.16 0 0.0 0
## 2 0 1.0 0 0.40 0 0.0 0
## 3 0 0.6 1 0.64 0 0.0 0
## 4 0 0.0 0 1.00 0 0.0 0
## 5 0 0.6 0 0.24 1 0.0 0
## 6 0 0.0 0 0.60 0 1.0 0
## 7 0 0.0 0 0.36 0 0.6 1

# Expected time to absorption
c = rep(1,7)
Ncmixed = Nmixed %*% c
Ncmixed

##   [,1]
## 1 1.56
## 2 1.40
## 3 2.24
## 4 1.00
## 5 1.84
## 6 1.60
## 7 1.96

# Absorption probabilities, B
Bmixed = Nmixed %*% Rmixed
Bmixed

##       0     8
## 1 0.936 0.064
## 2 0.840 0.160
## 3 0.744 0.256
## 4 0.600 0.400
## 5 0.504 0.496
## 6 0.360 0.640
## 7 0.216 0.784

# Probability of success, starting from 1 dollar:
p1mixed = Bmixed["1","8"]
p1mixed

## [1] 0.064

(b) Repeated matrix multiplications:

p=0.4
q=1-p
nodes=0:8
mixedmarkov

##     0 1   2 3   4 5   6 7   8
## 0 1.0 0 0.0 0 0.0 0 0.0 0 0.0
## 1 0.6 0 0.4 0 0.0 0 0.0 0 0.0
## 2 0.6 0 0.0 0 0.4 0 0.0 0 0.0
## 3 0.0 0 0.6 0 0.4 0 0.0 0 0.0
## 4 0.6 0 0.0 0 0.0 0 0.0 0 0.4
## 5 0.0 0 0.6 0 0.0 0 0.0 0 0.4
## 6 0.0 0 0.0 0 0.6 0 0.0 0 0.4
## 7 0.0 0 0.0 0 0.0 0 0.6 0 0.4
## 8 0.0 0 0.0 0 0.0 0 0.0 0 1.0

mixedresult = vector()
mixedtransient = vector()
epsilon = 1e-8
for (i in 1:10) {
  mixedpower = mixedmarkov %^% i

  # mixedmarkov raised to power i
  print(paste("i: ",i))
  print(round(mixedpower,8))
  
  mixedresult[i]=mixedpower["1","8"]
  mixedtransient[i] = sum(colSums(mixedpower)[2:7])
  if (i>1)
    mixeddiff=mixedresult[i]-mixedresult[i-1] 
  else mixeddiff=999

  print(paste(i,round(mixedtransient[i],10),round(mixedresult[i],   8),round(mixeddiff,10)))

  if (mixedtransient[i] < epsilon) {
    print(paste("Absorbed at power ", i))
    break
  }

}

## [1] "i:  1"
##     0 1   2 3   4 5   6 7   8
## 0 1.0 0 0.0 0 0.0 0 0.0 0 0.0
## 1 0.6 0 0.4 0 0.0 0 0.0 0 0.0
## 2 0.6 0 0.0 0 0.4 0 0.0 0 0.0
## 3 0.0 0 0.6 0 0.4 0 0.0 0 0.0
## 4 0.6 0 0.0 0 0.0 0 0.0 0 0.4
## 5 0.0 0 0.6 0 0.0 0 0.0 0 0.4
## 6 0.0 0 0.0 0 0.6 0 0.0 0 0.4
## 7 0.0 0 0.0 0 0.0 0 0.6 0 0.4
## 8 0.0 0 0.0 0 0.0 0 0.0 0 1.0
## [1] "1 3.6 0 999"
## [1] "i:  2"
##      0 1 2 3    4 5 6 7    8
## 0 1.00 0 0 0 0.00 0 0 0 0.00
## 1 0.84 0 0 0 0.16 0 0 0 0.00
## 2 0.84 0 0 0 0.00 0 0 0 0.16
## 3 0.60 0 0 0 0.24 0 0 0 0.16
## 4 0.60 0 0 0 0.00 0 0 0 0.40
## 5 0.36 0 0 0 0.24 0 0 0 0.40
## 6 0.36 0 0 0 0.00 0 0 0 0.64
## 7 0.00 0 0 0 0.36 0 0 0 0.64
## 8 0.00 0 0 0 0.00 0 0 0 1.00
## [1] "2 1 0 0"
## [1] "i:  3"
##       0 1 2 3 4 5 6 7     8
## 0 1.000 0 0 0 0 0 0 0 0.000
## 1 0.936 0 0 0 0 0 0 0 0.064
## 2 0.840 0 0 0 0 0 0 0 0.160
## 3 0.744 0 0 0 0 0 0 0 0.256
## 4 0.600 0 0 0 0 0 0 0 0.400
## 5 0.504 0 0 0 0 0 0 0 0.496
## 6 0.360 0 0 0 0 0 0 0 0.640
## 7 0.216 0 0 0 0 0 0 0 0.784
## 8 0.000 0 0 0 0 0 0 0 1.000
## [1] "3 0 0.064 0.064"
## [1] "Absorbed at power  3"

# Converged power matrix
mixedpower

##       0 1 2 3 4 5 6 7     8
## 0 1.000 0 0 0 0 0 0 0 0.000
## 1 0.936 0 0 0 0 0 0 0 0.064
## 2 0.840 0 0 0 0 0 0 0 0.160
## 3 0.744 0 0 0 0 0 0 0 0.256
## 4 0.600 0 0 0 0 0 0 0 0.400
## 5 0.504 0 0 0 0 0 0 0 0.496
## 6 0.360 0 0 0 0 0 0 0 0.640
## 7 0.216 0 0 0 0 0 0 0 0.784
## 8 0.000 0 0 0 0 0 0 0 1.000

# Sums of columns
print("Column sums: ")

## [1] "Column sums: "

print(round(colSums(mixedpower),8))

##   0   1   2   3   4   5   6   7   8 
## 5.2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.8

# mixed result
print(paste("mixed result: ",round(mixedresult[i],8)))

## [1] "mixed result:  0.064"

For this strategy, the probability of winning is 0.064 .

How does this probability compare with that obtained for the bold strategy?

The result is the same - for both strategies, the probability of success is 6.4% .

The only difference between the two strategies occurs if the prisoner starts with three dollars:

In the bold strategy, you would bet 3 dollars when you have 3 dollars
In the mixed strategy, you would bet 1 dollar when you have 3 dollars

In all other cases, the betting is unchanged.

It is not possible to enter state 3 unless you start there, so for this specific case, it becomes moot.
But, even if you did start from that state, the probability results of both strategies are the same - no matter which state you start from.

There can be no more than 3 plays, no matter which state you start from.

CH11_S2_EX13b14_P424

Michael Y.

10/30/2019