#
# Two first order Markov chain with 5 states
#
# P1 - transition matrix with good mixing
# P2 - transition matrix with bad mixing

P1 = matrix(c(
0.35,0.35,0.10,0.10,0.10,
0.15,0.55,0.10,0.10,0.10,
0.15,0.15,0.10,0.20,0.40,
0.15,0.15,0.40,0.10,0.20,
0.15,0.15,0.20,0.40,0.10),5,5,byrow=TRUE)

eps = 0.001

P2 = matrix(c(
0.4985,0.4985,0.001,0.001,0.001,
0.2970,0.7000,0.001,0.001,0.001,
0.0015,0.0015,0.199,0.299,0.499,
0.0015,0.0015,0.499,0.199,0.299,
0.0015,0.0015,0.299,0.499,0.199),5,5,byrow=TRUE)

# P2 = matrix(c(
# 0.499850,0.499850,0.000100,0.000100,0.000100,
# 0.299850,0.699850,0.000100,0.000100,0.000100,
# 0.000150,0.000150,0.199900,0.299900,0.499900,
# 0.000150,0.000150,0.499900,0.199900,0.299900,
# 0.000150,0.000150,0.299900,0.499900,0.199900),5,5,byrow=TRUE)

prod = function(P,k){
  PP = P
  for (i in 1:k)
    PP = PP%*%P
  return(PP) 
}

prod(P1,2^5)
prod(P2,2^5)
prod(P2,2^10)
prod(P2,2^16)

# Equilibrium distribution
pi1.true = prod(P1,2^5)[1,]
pi2.true = prod(P2,2^16)[1,]

# Simulaton the behaviors of both Markov chains
M0    = 200000
M     = 20000
#M0    = 1000
#M     = 1000
niter = M0+M
P     = P1
xs    = rep(0,niter)
xs[1] = 1
for (i in 2:niter)
 xs[i] = sample(1:5,prob=P[xs[i-1],],size=1)
xs1 = xs[(M0+1):niter]
P     = P2
xs    = rep(0,niter)
xs[1] = 5
for (i in 2:niter)
 xs[i] = sample(1:5,prob=P[xs[i-1],],size=1)
xs2 = xs[(M0+1):niter]

pi1 = table(xs1)/M
pi2 = table(xs2)/M

par(mfrow=c(2,2),mai=c(0.8,0.8,0.3,0.3))
plot(xs1,ylim=c(1,5),xlab="Iterations of the Markov chain",ylab="States")
title("Transition matrix: P1")
plot(xs2,ylim=c(1,5),xlab="Iterations of the Markov chain",ylab="States")
title("Transition matrix: P2")
plot(pi1,type="h",ylim=range(pi,pi1,pi2),xlab="States",ylab="Equilibrium distribution",lwd=3)
lines(1:5+0.1,pi1.true,type="h",col=2,lwd=3)
plot(pi2,type="h",ylim=range(pi,pi1,pi2),xlab="States",ylab="Equilibrium distribution",lwd=3)
lines(1:5+0.1,pi2.true,type="h",col=2,lwd=3)
legend("topright",legend=c("True","Estimated"),col=2:1,lwd=3,bty="n")