# 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')