### Parallel Tempering

We want to do inference on a target density $$\pi$$ and make use of parallel cores. We also want to escape from local modes in cases where there are a few. For instance introduce the following toy target density: $\pi(x) \propto \exp\left(- 10 (x^2 - 1)^2\right).$ The target density has a mode at $$1$$ and a mode at $$-1$$. We introduce the tempered density $\pi(x)^\gamma \propto \exp\left(- 10 \gamma (x^2 - 1)^2\right)$ for any $$\gamma\in [0,1]$$. Let us look at the target density function for $$\gamma = 0.1, 0.4, 0.7, 1$$.

target <- function(x) -10 * (x * x - 1) * (x * x - 1)

for (temp in seq(from = 0.1, to = 1, length.out = 4)) {
cst <- integrate(f = function(x) exp(temp * target(x)), lower = -4, upper = 3)\$value
norm_target <- function(x) exp(temp * target(x))/cst
g <- qplot(x = c(-5, 5), geom = "blank") + stat_function(fun = function(x) norm_target(x)^(temp),
colour = "red", size = 1, n = 1000) + xlim(-3, 3) + xlab(paste0("gamma = ",
temp)) + ylab("density") + ylim(0, 2)
print(g)
}

We introduce $$N$$ values of $$\gamma$$ between $$0.1$$ and $$1$$. We can run an independent MCMC for each $$\gamma$$.

niterations <- 10000
nchains <- 10
inv_temperatures <- seq(from = 0.1, to = 1, length.out = nchains)
current_values <- rnorm(nchains, sd = 10)
current_pdf <- inv_temperatures * target(current_values)
chains <- matrix(ncol = nchains, nrow = niterations)
chains[1, ] <- current_values
for (t in 2:niterations) {
proposals_values <- rnorm(nchains, mean = current_values, sd = 0.5)
proposals_pdf <- inv_temperatures * target(proposals_values)
accepts <- (log(runif(nchains)) < (proposals_pdf - current_pdf))
current_pdf[accepts] <- proposals_pdf[accepts]
current_values[accepts] <- proposals_values[accepts]
chains[t, ] <- current_values
}
df <- melt(chains)

For the chains corresponding to high values of $$\gamma$$, we see that the chains struggle to move from one mode to the other.

gtraceplot1 <- qplot(data = subset(df, X2 >= 8), x = X1, y = value, colour = factor(X2),
geom = "line") + ylab("X") + xlab("iteration") + ylim(-2, 2)
gtraceplot1 <- gtraceplot1 + theme(legend.position = "bottom")
gtraceplot1 <- gtraceplot1 + scale_color_discrete(name = "# chain")
print(gtraceplot1)