#------------------------------------------------------ # # Metropolis algorithm on three states # # GdA 4/4/2016 (rev. 10/5/2010) #------------------------------------------------------ p <- c(0.2, 0.3, 0.5) N <- length(p) n = 10000 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[1:100], pch=19, ty='b', col='blue', cex=0.4) par(mar=old.mar)