Write a program to compute \(u^{(n)}\) given u and P. Use this program tocompute \(u^{(10)}\) for the Land of Oz example, with \(u = (0, 1, 0)\), and with \(u = (1/3, 1/3, 1/3)\).

Note: see page probability text ex 11.1 (p. 406) for the details of “The Land of Oz” as referred to in the question.

#the function as per the question
library(expm)

compute.un <- function(n,u,P){
  
  #could be done in 1 line... kept it separate for my own benefit
  
  Pn <- P %^% n #compute p^n
  
  out <- u %*% Pn  #compute uP^n

  return(out)
}
#set up the transition max for oz
oz <- matrix(c(0.5,0.25,0.25,0.5,0,0.5,0.25,0.25,0.5),nrow = 3,byrow=T)

u1 <- c(0,1,0)
u2 <- c(0.33,0.33,0.33)

#for u1
round(compute.un(10,u1,oz),2)
##      [,1] [,2] [,3]
## [1,]  0.4  0.2  0.4
#for u2
round(compute.un(10,u2,oz),2)
##      [,1] [,2] [,3]
## [1,]  0.4  0.2  0.4

Despite the fact that I have the same output for both, I’m cautiously optimistic that this is correct.