#------------------------------------------------------
#
#  Metropolis algorithm on three states
#
#  GdA 4/4/2016
#------------------------------------------------------

 p <- c(0.2, 0.3, 0.5)
N <- length(p)

n = 200
x <- rep(0, n)
x[1] = 1

for (i in 2:n) {
   x.p <- sample(1:N, 1)
   A   <- min( 1, p[x.p]/p[x[i-1]] )
   x[i] = ifelse ( runif(1) <= A, x.p, x[i-1] )
}

ns = 10
print(table(x[(ns+1):n])/(n-ns))   # skip the first ns steps

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)