2

This is my task:

Peter goes to the Casino with 1 dollar. With the chance of p, Peter wins 1 dollar and the chance of (1-p) he looses 1 dollar. The process can be seen as a markov chain.

If Peter reaches 0 dollars he goes home bankrupt, if he manages to reach 5 dollars he goes home happy.

Find the probability of Peter going home with 5 dollars when p= 30%, 40%, 50%, 60% & 70%. Construct matrixes for each probability where the first 4 states are the transient class( 1- 4 dollars) and the last two states are the two recurrent states( 0 & 5 dollars).

My plan in solving it

Find when each separate matrix converges(P^n = P^n+1) with when_converged.

Then use that n in mpow to see the probability of going from 1 dollars to 5 dollars, in other words from state 1 to 6.

This is my code:

mpow <- function(P, n) {
    if (n == 0) {
        return(diag(nrow(P)))
    } else if (n == 1) {
        return(P)
    } else {
        return(P %*% mpow(P, n - 1))
    }
}

when_converged <- function(P, tol=0.00005) { 
    n = 1; diff = 1 
    while (diff > tol) { 
        A <- mpow(P, n) 
        B <- mpow(P, n+1) 
        diff <- max(abs(A - B)) 
        n <- n + 1 
    } 
    return(n) 
}


P30 <- matrix(c(0, 0.3, 0, 0, 0.7, 0, 0.7, 0, 0.3, 0, 0, 0, 0, 0.7, 0, 0.3, 0, 0, 0, 0, 0.7, 0, 0, 0.3, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1), nrow = 6, ncol = 6, byrow = TRUE)

P40 <- matrix(c(0, 0.4, 0, 0, 0.6, 0, 0.6, 0, 0.4, 0, 0, 0, 0, 0.6, 0, 0.4, 0, 0, 0, 0, 0.6, 0, 0, 0.4, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1), nrow = 6, ncol = 6, byrow = TRUE)

P50 <- matrix(c(0, 0.5, 0, 0, 0.5, 0, 0.5, 0, 0.5, 0, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1), nrow = 6, ncol = 6, byrow = TRUE)

P60 <- matrix(c(0, 0.6, 0, 0, 0.4, 0, 0.6, 0, 0.4, 0, 0, 0, 0, 0.6, 0, 0.4, 0, 0, 0, 0, 0.6, 0, 0, 0.4, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1), nrow = 6, ncol = 6, byrow = TRUE)

P70 <- matrix(c(0, 0.7, 0, 0, 0.3, 0, 0.7, 0, 0.3, 0, 0, 0, 0, 0.7, 0, 0.3, 0, 0, 0, 0, 0.7, 0, 0, 0.3, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1), nrow = 6, ncol = 6, byrow = TRUE)

when_converged(P30, 0.00005)

From Rstudio I get that P30 is converged at 35.

when_converged(P40, 0.00005)

From Rstudio I get that P40 is converged at 37.

when_converged(P50, 0.00005)

From Rstudio I get that P50 is converged at 47.

when_converged(P60, 0.00005)

From Rstudio I get that P60 is converged at 61.

when_converged(P70, 0.00005)

From Rstudio I get that P70 is converged at 79.

mpow(P30, 35)

mpow(P40, 37)

mpow(P50, 47)

mpow(P60, 61)

mpow(P70, 79)

What I need help with

What I get from Rstudio is that for mpow(P60, 61) & mpow(P70, 79) the probability of going home with 5 dollars becomes less compared to mpow(P50, 47) & mpow(P40, 37). Where the probability of winning 1 dollar is less. Which feels wrong. Is there anything I'm doing wrong? Please try to solve it using my method & not with an entire different code.

IRTFM
  • 258,963
  • 21
  • 364
  • 487
PeterNiklas
  • 75
  • 1
  • 9
  • If anything is unclear, please ask. – PeterNiklas Jul 11 '16 at 17:57
  • 1
    As I see it the probability of remaining in either the 0 or the 5 state should be 1. That's what is called an absorbing state. There should only be transitions from a state of one to either 0 or 2 and they should be complementary (i.e. sum to one). Your matrices do no look like that to me. – IRTFM Jul 12 '16 at 00:14

1 Answers1

3

This is how I constructed a P30 matrix ... not the same as yours:

> P30 <- matrix(c(1,  0, 0,  0,  0, 0,
+                   0.7, 0, 0.3, 0, 0,  0, 
+                   0,  0.7, 0, 0.3, 0, 0, 
+                   0,  0, 0.7, 0, 0.3, 0,
+                   0,  0,  0,  0.7, 0 , 0.3,                  
+                   0,  0, 0, 0, 0, 1), nrow = 6, ncol = 6, byrow = TRUE)
> P30
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]  1.0  0.0  0.0  0.0  0.0  0.0
[2,]  0.7  0.0  0.3  0.0  0.0  0.0
[3,]  0.0  0.7  0.0  0.3  0.0  0.0
[4,]  0.0  0.0  0.7  0.0  0.3  0.0
[5,]  0.0  0.0  0.0  0.7  0.0  0.3
[6,]  0.0  0.0  0.0  0.0  0.0  1.0

Notice that in each row the input column gets sent only to itself in the case of the 0 or 5 state, but in the others it gets sent to an adjacent output column. So in-1 goes to either out-0 or out-2. Probably clearer to display with col- and row-names:

> rownames(P30) <- 0:5
> colnames(P30) <- 0:5
> P30
    0   1   2   3   4   5
0 1.0 0.0 0.0 0.0 0.0 0.0
1 0.7 0.0 0.3 0.0 0.0 0.0
2 0.0 0.7 0.0 0.3 0.0 0.0
3 0.0 0.0 0.7 0.0 0.3 0.0
4 0.0 0.0 0.0 0.7 0.0 0.3
5 0.0 0.0 0.0 0.0 0.0 1.0

This could aid in building such matrices with different values for P

 p0 <-  matrix(0, nrow = 6, ncol = 6); p=.30
 p30 <- p0; p30 [cbind(2:5,1:4)] <- 1-p
             p30[cbind(2:5,3:6)] <- p
 p30[ cbind(c(1,1),c(6,6))] <- 1
 p30

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]  1.0  0.0  0.0  0.0  0.0  0.0
[2,]  0.7  0.0  0.3  0.0  0.0  0.0
[3,]  0.0  0.7  0.0  0.3  0.0  0.0
[4,]  0.0  0.0  0.7  0.0  0.3  0.0
[5,]  0.0  0.0  0.0  0.7  0.0  0.3
[6,]  0.0  0.0  0.0  0.0  0.0  1.0

The deterministic or theoretical probabilities after three iterations (starting at state=1:

c(0,1,0,0,0,0) %*% P30 %*% P30 %*% P30
#-----
         0 1     2 3     4 5
[1,] 0.847 0 0.126 0 0.027 0

Agrees with your mpow

> c(0,1,0,0,0,0) %*% mpow(P30 ,3)
         0 1     2 3     4 5
[1,] 0.847 0 0.126 0 0.027 0

There's also a matrix-power function %^% in the expm-package.

> c(0,1,0,0,0,0) %*% expm::'%^%'( P30,3)
         0 1     2 3     4 5
[1,] 0.847 0 0.126 0 0.027 0
IRTFM
  • 258,963
  • 21
  • 364
  • 487