# file "Rtempered" -- a tempered Metropolis algorithm in one dimension

M = 101000  # run length
B = 1000  # amount of burn-in
maxtemp = 10
temp = 1  # 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,4); plot(tf,-10,30)
# tf = function(x) pi(x,8); 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")