# file "Rtempered" -- a tempered Metropolis algorithm in one dimension M = 51000 # run length B = 1000 # amount of burn-in maxtemp = 10 temp = 4 # initial value, if no tempering dotempering = TRUE pi = function(x, thetemp) { if ( (thetemp<1) || (thetemp>maxtemp) ) return(0.0) else return( 0.5 * dnorm(x, 0, thetemp) + 0.5 * dnorm(x, 20, thetemp) ); } h = function(y) { return(y) } xlist = rep(0,M) # for keeping track of chain values templist = rep(0,M) # for keeping track of temperature values hlist = rep(0,M) # for keeping track of h function values numxaccept = numtempaccept = temptries = 0; sigma = 1 # proposal scaling X = runif(1,-10,30) # overdispersed starting distribution if (dotempering) temp = floor(1+runif(1,0,maxtemp)) # overdispersed starting distribution for (i in 1:M) { # PROPOSED X MOVE Y = X + sigma * rnorm(1) # proposal value U = runif(1) # for accept/reject A = pi(Y,temp) / pi(X,temp) # for accept/reject if (U < A) { X = Y # accept proposal numxaccept = numxaccept + 1; } # PROPOSED TEMP MOVE if (dotempering) { temptries = temptries + 1 newtemp = temp - 1 + 2 * floor( runif(1,0,2) ) # temp +- 1 U = runif(1) # for accept/reject A = pi(X,newtemp) / pi(X,temp) # for accept/reject if (U < A) { temp = newtemp # accept proposal numtempaccept = numtempaccept + 1; } } xlist[i] = X; templist[i] = temp; hlist[i] = h(X); } countedlist = c( rep(FALSE,B), rep(TRUE,M-B) ) valslist = hlist[(templist==1) & (countedlist==TRUE)] allvalslist = hlist[(countedlist==TRUE)] cat("Ran tempered algorithm for", M, "iterations, with burn-in", B, "\n"); cat("x acceptance rate =", numxaccept/M, "\n"); cat("temp acceptance rate =", numtempaccept/M, "\n"); cat("number of sampled values =", length(valslist), "\n") cat("mean of sampled values is about", mean(valslist), "\n") if (length(valslist) > 1) { se1 = sd(valslist) / sqrt(length(valslist)) cat("iid standard error would be about", se1, "\n") varfact <- function(xxx) { 2 * sum(acf(xxx, plot=FALSE)$acf) - 1 } cat("true standard error is about", se1 * sqrt( varfact(valslist) ), "\n") } cat("[ mean of ALL values is about", mean(allvalslist), "]\n") # tf = function(x) pi(x,1); plot(tf,-10,30) # tf = function(x) pi(x,10); plot(tf,-10,30) plot(xlist, type='l') # points((1:M)[templist==1], xlist[templist==1], col="red")