## Calculating n-step transition probability matrices: ######################################################## # function for calculating matrix powers, by Alberto Monteiro # Copy and paste this into R first: # matpowfast <- function(mat, n) { if (n == 1) return(mat) result <- diag(1, ncol(mat)) while (n > 0) { if (n %% 2 != 0) { result <- result %*% mat n <- n - 1 } mat <- mat %*% mat n <- n / 2 } return(result) } # ######################################################## # Example 1, Section 4.2 notes: P = matrix(c(0.5,0.5,0.3,0.7), nrow=2, ncol=2, byrow=T) P.5 <- matpowfast(P,5) print(P.5) # Example 3, Section 4.2 notes: P = matrix(c(0,1,0,0.25,0.5,0.25,0,1,0), nrow=3, ncol=3, byrow=T) P.4 <- matpowfast(P,4) print(P.4) # Example 4, Section 4.2 notes: Q = matrix(c(7/12,5/12,0,0, 7/12,0,5/12,0, 7/12,0,0,5/12, 0,0,0,1), nrow=4, ncol=4, byrow=T) Q.10 <- matpowfast(Q,10) print(Q.10) Q.9 <- matpowfast(Q,9) print(Q.9) # Example 1, limiting probabilities: P = matrix(c(0.5,0.5,0.3,0.7), nrow=2, ncol=2, byrow=T) print(matpowfast(P,4)) print(matpowfast(P,8)) print(matpowfast(P,12)) ######################################################## # Roulette Examples, Section 4.6 notes: # Setting up the P.T matrix: P.T <- matrix(0,nrow=9,ncol=9) for (i in 1:8){ P.T[i,i+1] <- 18/38 P.T[i+1,i] <- 20/38 } I9 <- diag(rep(1,times=9)) solve(I9-P.T) sum(solve(I9-P.T)[6,]) ######################################################## # Extinction probability, Section 4.7 notes: # Putting all the terms on one side of the equation, we have: # 0 = 0.3 - 0.555*pi0 + 0.165*pi0^2 + 0.060*pi0^3 + 0.022*pi0^4 + 0.008*pi0^5, so: coef.vector <- c(0.3, -0.555, 0.165, 0.060, 0.022, 0.008) # coefficients must be listed in order: # first the constant term, then the lowest-degree term to the highest-degree term. polyroot(coef.vector) # The smallest positive real root is 0.841088.