---
title: "Double well"
author: "Bob Verity and Pete Winskill"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Double well}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, echo = FALSE}
# set random seed
set.seed(1)

# load the drjacoby package
library(drjacoby)
```

**Purpose:** to compare *drjacoby* results for a challenging problem involving a multimodal posterior, both with and without temperature rungs.

## Model

We assume a single parameter `mu` drawn from a double well potential distribution, defined by the formula:

$$
\begin{aligned}
  \mu &\propto exp\left(-\gamma(\mu^2 - 1)^2\right)
\end{aligned}
$$
where $\gamma$ is a parameter that defines the strength of the well (higher $\gamma$ leads to a deeper valley and hence more challenging problem). NB, there is no data in this example, as the likelihood is defined exactly by these parameters.

Likelihood and prior:

```{r, echo = FALSE, comment = ''}
Rcpp::sourceCpp(system.file("extdata/checks/", "doublewell_loglike_logprior.cpp", package = 'drjacoby', mustWork = TRUE))
```

Parameters dataframe:

```{r}
L <- 2
gamma <- 30
df_params <- define_params(name = "mu", min = -L, max = L,
                           name = "gamma", min = gamma, max = gamma)
```


## Single temperature rung (no Metropolis coupling)

```{r}
mcmc <- run_mcmc(data = list(x = -1),
                 df_params = df_params,
                 loglike = "loglike",
                 logprior = "logprior",
                 burnin = 1e3,
                 samples = 1e5,
                 chains = 1,
                 rungs = 1,
                 silent = TRUE)
```

```{r}
# trace plot
plot_trace(mcmc, show = "mu")
```

```{r}
# extract posterior draws
output_sub <- subset(mcmc$output, phase == "sampling")
mu_draws <- output_sub$mu

# get analytical solution
x <- seq(-L, L, l = 1001)
fx <- exp(-gamma*(x^2 - 1)^2)
fx <- fx / sum(fx) * 1/(x[2]-x[1])

# overlay plots
hist(mu_draws, breaks = seq(-L, L, l = 201), probability = TRUE, main = "", col = "black")
lines(x, fx, col = 2, lwd = 2)
```


## Multiple temperature rungs

```{r}
mcmc <- run_mcmc(data = list(x = -1),
                 df_params = df_params,
                 loglike = "loglike",
                 logprior = "logprior",
                 burnin = 1e3,
                 samples = 1e5,
                 chains = 1,
                 rungs = 11,
                 alpha = 2,
                 pb_markdown = TRUE)
```

```{r, fig.width=6, fig.height=5}
# trace plot
plot_trace(mcmc, show = "mu")

# coupling acceptance plot
plot_mc_acceptance(mcmc)
```

```{r, fig.width=6, fig.height=5}
# extract posterior draws
output_sub <- subset(mcmc$output, phase == "sampling")
mu_draws <- output_sub$mu

# overlay plots
hist(mu_draws, breaks = seq(-L, L, l = 201), probability = TRUE, main = "", col = "black")
lines(x, fx, col = 2, lwd = 2)
```