# hhh vs. thh game setwd("/Users/pollard/Academic/DP-Courses/251/251.spring2013/Handouts/") states <- c("nix","H","T","HH", "TH","THH","HHH") # create the transition matrix P <- matrix(0,nrow=7,ncol=7) # first fill with zeros rownames(P) <- colnames(P) <- states P["nix","H"] <- P["nix","T"] <- P["H","T"] <- P["H","HH"] <- P["T","T"] <- 1/2 P["T","TH"] <- P["TH","THH"] <- P["TH","T"] <- P["HH","HHH"] <- P["HH","T"] <- 1/2 P["HHH","HHH"] <- P["THH","THH"] <-1 # Solve theta = P %*% theta, with theta["HHH"]=1 and theta["THH"]=0 I7 <- diag(7) # the 7x7 identity matrix constraints <- I7[6:7,] M <- I7 - P # Want to solve M %*% theta = 0 subject to constraints %*% theta = c(0,1) # Note that the HHH and THH rows are now useless: filled with zeros. # Replace them by the constraints fixrows <- c("THH","HHH") M[fixrows, ] <- constraints theta <- solve(M,c(0,0,0,0,0,0,1)) names(theta) <- states # probs of reaching "HHH" for various starting states # Compare with: ####################################################### mpower = function(mm,kk){#raise matrix mm to power 2^kk for (i in 1:kk){ mm = mm %*% mm } return(mm) } ####################################################### round(mpower(P,10),3) # raise P to the power of 2^10 # Expected duration tau <- solve(M,c(1,1,1,1,1,0,0)) names(tau) <- states # Solutions: # > theta # nix H T HH TH THH HHH # 0.125 0.250 0.000 0.500 0.000 0.000 1.000 # > tau # nix H T HH TH THH HHH # 7 6 6 4 4 0 0