# file "Rpara" -- a parallel-tempered Metropolis algorithm in one dimension M = 10100 # run length B = 100 # amount of burn-in maxtemp = 8 dotempering = TRUE pi = function(x) { return( 0.5 * dnorm(x, 0, 1) + 0.5 * dnorm(x, 20, 1) ); } g = function(x, thetemp) { if ( (thetemp<1) || (thetemp>maxtemp) ) return(0.0) else return( pi(x)^(1/thetemp) ) } h = function(y) { return(y) } xlist = matrix( rep(0,M*maxtemp), ncol=maxtemp ) # for chain values numxaccept = rep(0,maxtemp) numtempaccept = 0; sigma = 5 # proposal scaling X = runif(maxtemp,-10,30) # overdispersed starting distribution, dim=maxtemp for (i in 1:M) { for (temp in 1:maxtemp) { # PROPOSED X[temp] MOVE Y = X[temp] + sigma * rnorm(1) # proposal value U = runif(1) # for accept/reject A = g(Y,temp) / g(X[temp],temp) # for accept/reject if (U < A) { X[temp] = Y # accept proposal numxaccept[temp] = numxaccept[temp] + 1; } } # PROPOSED TEMP SWAP if (dotempering) { j = floor(1+runif(1,0,maxtemp)) # uniform on {1,2,...,maxtemp} k = floor(1+runif(1,0,maxtemp)) # uniform on {1,2,...,maxtemp} # j = floor(1+runif(1,0,maxtemp-1)) # uniform on {1,2,...,maxtemp-1} # k = j + 1; U = runif(1) # for accept/reject A = (g(X[j],k) * g(X[k],j)) / (g(X[j],j) * g(X[k],k)); if (U < A) { # accept proposed swap tmpval = X[j]; X[j] = X[k]; X[k] = tmpval; numtempaccept = numtempaccept + 1; } } xlist[i,] = X; } valslist = h(xlist[((B+1):M),1]) cat("Ran parallel 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("mean of sampled values is about", mean(valslist), "\n") 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") plot(xlist[,1], type='l') # plot(xlist[(1000:2000),1], type='l') # plot(xlist[(1000:2000),4], type='l') # plot(xlist[(1000:2000),8], type='l')