updating, more to come

Author: maxbiostat

assignments/A1.pdf

@@ -167,7 +167,7 @@ chains_coupled = function(pot = U,
   xmat
 }
 
-mc3 <- chains_coupled(tune = 1)
+mc3 <- chains_coupled(tune = .1)
 
 mcmcSummary(mc3, rows = length(temps))
 
@@ -185,4 +185,4 @@ for (k in 1:ncol(mc3)){
 mean(theta.draws) ## MCMC (HMC)
 posterior.mean.quadrature ## quadrature
 m_is ## Importance sampling
-colMeans(mc3)
+colMeans(mc3)

assignments/A1_solutions.pdf

@@ -45,87 +45,3 @@ curve(
   main = "Cauchy example -- BC exercise 1.28"
 )
 
-############ 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)
-}
-