#------------------------------------------------------
#
#  Three state MCMC (cfr Andrieu at al.)
#
#  GdA 2/4/2016
#------------------------------------------------------

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

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

T <- rbind( c(0, 1, 0), c(0, 0.1, 0.9), c(0.6, 0.4, 0) )
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
print(table(x[(ns+1):n])/(n-ns))   # skip the first ns steps

# 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', cex=0.4)
par(mar=old.mar)


cat("\n")
print(x[1:30])