# Markov Chain code

pp<-matrix(c(.1,.2,.7,.2,.2,.6,.6,.1,.3),ncol=3,byrow=T)

ipp<-diag(1,3)-t(pp)
tipp<-ipp
tipp[2,]<-rep(1,3)
ypp<-c(0,1,0)
solve(tipp,ypp)

z<-numeric(100000)
z[1]<-1
for(i in 2:100000){
	z[i] <- sample(c(1,2,3),prob=pp[z[i-1],],size=1)
}

z.converged <- z[11:100000]
length(z.converged[z.converged==1])/length(z.converged)
# or mean(z.converged==1)
length(z.converged[z.converged==2])/length(z.converged)
# or mean(z.converged==2)
length(z.converged[z.converged==3])/length(z.converged)
# or mean(z.converged==3)

hist(z.converged)

pp2 <- pp%*%pp
pp4 <- pp2%*%pp2
pp8 <- pp4%*%pp4
pp16 <- pp8%*%pp8
pp32 <- pp16%*%pp16
pp64 <- pp32%*%pp32

pp64