all in one place

Author: maxbiostat

code/cauchy_bimodal.r

@@ -97,7 +97,92 @@ hist(theta.draws, probability = TRUE,
 curve(c_posterior, lwd = 3, add = TRUE)
 lines(density(IS.draws), lwd = 3, lty = 2, col = 2)
 
+
+############ Parallel tempering
+
+U = function(gam, x)
+{
+  - gam * c_logposterior_kernel(x)
+}
+
+curried = function(gam)
+{
+  message(paste("Returning a function for gamma =", gam))
+  function(x)
+    U(gam, x)
+}
+U4 = curried(4)
+
+op = par(mfrow = c(2, 1))
+curve(U4(x), minT, maxT, main = "Potential function, U(x)")
+curve(exp(-U4(x)), minT, maxT, main = "Unnormalised density function, exp(-U(x))")
+par(op)
+
+
+chain = function(target, tune = 0.1, init = 1)
+{
+  x = init
+  xvec = numeric(iters)
+  for (i in 1:iters) {
+    can = x + rnorm(1, 0, tune)
+    logA = target(x) - target(can)
+    if (log(runif(1)) < logA)
+      x = can
+    xvec[i] = x
+  }
+  xvec
+}
+
+temps = 2 ^ (0:3)
+iters = 1e5
+
+mat = sapply(lapply(temps, curried), chain)
+colnames(mat) = paste("gamma=", temps, sep = "")
+
+require(smfsb)
+mcmcSummary(mat, rows = length(temps))
+
+
+chains_coupled = function(pot = U,
+                          tune = 0.1,
+                          init = 1)
+{
+  x = rep(init, length(temps))
+  xmat = matrix(0, iters, length(temps))
+  for (i in 1:iters) {
+    can = x + rnorm(length(temps), 0, tune)
+    logA = unlist(Map(pot, temps, x)) - unlist(Map(pot, temps, can))
+    accept = (log(runif(length(temps))) < logA)
+    x[accept] = can[accept]
+    # now the coupling update
+    swap = sample(1:length(temps), 2)
+    logA = pot(temps[swap[1]], x[swap[1]]) + pot(temps[swap[2]], x[swap[2]]) -
+      pot(temps[swap[1]], x[swap[2]]) - pot(temps[swap[2]], x[swap[1]])
+    if (log(runif(1)) < logA)
+      x[swap] = rev(x[swap])
+    # end of the coupling update
+    xmat[i, ] = x
+  }
+  colnames(xmat) = paste("gamma=", temps, sep = "")
+  xmat
+}
+
+mc3 <- chains_coupled(tune = 1)
+
+mcmcSummary(mc3, rows = length(temps))
+
+
+par(mfrow = c(2, 2))
+for (k in 1:ncol(mc3)){
+  hist(mc3[, k], probability = TRUE,
+       xlim = c(minT, maxT),
+       main = paste("gamma =", temps[k]),
+       xlab = expression(theta))
+  curve(c_posterior, lwd = 3, add = TRUE)
+}
+
 ##### Posterior mean estimates
 mean(theta.draws) ## MCMC (HMC)
 posterior.mean.quadrature ## quadrature
-m_is ## Importance sampling
+m_is ## Importance sampling
+colMeans(mc3)