#------------------------------------------------------
#
#  Three state MCMC (cfr Andrieu at al.)
#
#  GdA 2/4/2016, rev. 10/5/2020
#------------------------------------------------------

extract.cp <- function(cp) {  # extract using cumulative distribution
  r <- runif(1)
  for (i in 1:3) {
    if(r <= cp[i]) return(i)
  }
} 

n = 200  # 200
x <- rep(0, n)

# transition matrix
T <- rbind( c(0, 1, 0), c(0, 0.1, 0.9), c(0.6, 0.4, 0) )
# cumulative prob. distr. (-> for random -> extract.cp)
cp <- rbind( cumsum(T[1,]),   cumsum(T[2,]), cumsum(T[3,]) )


x[1] = 1

for (i in 2:n) {
   x[i] <- extract.cp(cp[x[i-1],])
}

ns = 10  # number of steps to be skipped (to lose memory)

print(table(x[(ns+1):n]))          # freq. of visiting the 3 states
print(table(x[(ns+1):n])/(n-ns))   # relative frequency 

# commands to perform partial checks (n=10100)
# y=x[101:10100] 
# for (i in 1:10) { print(table(y[((i-1)*1000+1):(i*1000)])/1000); print(table(y[1:(i*1000)])/(i*1000))} 

old.mar = par("mar")
par(mar=c(4,4,0.1,0.1))
plot(x, pch=19, col='blue', ty='b', lty=1, lwd=0.5, cex=0.8)
par(mar=old.mar)


cat("\n First 37 moves\n")
print(x[1:37])