---
title: 'Module 9: Solutions'
author: ""
date: ""
output: html_document
---

## Exercise 1: Minimise auto-correlation

Write a function in R that uses a Metropolis-Hastings algorithm to simulate from
a standard normal distribution, using normally distributed proposals centred at
the current value (in other words, this is a random walk Metropolis algorithm).
Make sure your code allows you to adjust the standard deviation of the proposal
distribution and monitor the acceptance probability. (**Hint:** You already
implemented this in the exercises for module 7 (see the solution if you didn't
finish the exercise) and you may use this as the main part of your function).
The skeleton of such a function is given below:
```{r}
library(coda)
myMH <- function(N, sigma){
  x <- rep(0, N) # Empty Markov chain
  a <- 0 # Number of accepted proposals
  # Below implement the M-H alg. and add 1 to `a` every time you accept
  x[1] <- 0 ## initial value
  for(i in 2:N){
    y <- rnorm(1, x[i-1], sigma)
    u <- runif(1)
    H <- dnorm(y)/dnorm(x[i-1])
    if(u<H){
      x[i] <- y
      a <- a+1
    }else{
      x[i] <- x[i-1]
    }
  }
  # Finally return a list of the results
  return(list(x=mcmc(x), a=a/(N-1)))
}
# A short run with proposal std. dev. 1
rslt <- myMH(10000, 1)
```

To plot the auto-correlation function use the function `acf()` in
R. To approximately calculate $V = 1 + 2\sum_{m=1}^\infty \rho_m$ use the
following command in R:
```{r}
V <- 2 * sum(acf(rslt$x, plot = FALSE)$acf) - 1
```

Make some experiments and find the standard deviation for the proposal distribution that gives
the smallest value of $V$. What is the acceptance probability
in this case? (Make sure that you have run the MH-algorithm long
enough (really long), so that the auto-correlation function is well estimated.)
```{r}
sigma <- seq(1, 5, by = .5)
V <- rep(0, length(sigma))
a <- rep(0, length(sigma))
for(i in 1:length(sigma)){
  rslt <- myMH(100000, sigma[i])
  x <- rslt$x
  acf_val <- acf(x, plot = FALSE)$acf
  V[i] <- 2 * sum(acf_val) - 1
  a[i] <- rslt$a
}
plot(sigma, V, type = "l")
a[which.min(V)]
```


## Exercise 2: Beetles

This example is concerned with different beetles exposed in 5 hours to different
concentrations of carbon disulphide. The data contains for 481 beetles, in the
first column, the concentration, and in the second column, the status of the
beetles after the exposure. The status is coded as zeros (dead beetles) and ones 
(survived beetles). Read in the data set `beetles.dat` from the course website:
```{r}
data_url <- "https://asta.math.aau.dk/course/bayes/2021/?file=beetles.dat"
beetles <- read.table(data_url)
```

Let $x_i>0$ and $y_i\in\{0,1\}$ denote the concentration and status for the $i$th
beetle, respectively. We use a logistic regression model for the data of survived and dead beetles, i.e., we condition on the concentrations and assume
that the probability that the $i$th beetle survives its given dose $x_i$ is
$$
\frac{\exp(a+bx_i)}{1+\exp(a+bx_i)}
$$
where $a$ and $b$ are real parameters.
In R you can define the data vecors and log-likelihood function as:
```{r}
beet_x <- beetles[,1]
beet_y <- beetles[,2]
llik = function(a, b){
  ## BEWARE: We refer to global variables beet_x and beet_y here!
  return(sum(beet_y*(a+b*beet_x)) - sum(log(1+exp(a+b*beet_x))))
} 
```


Now, let $\text{Cauchy}(\gamma)$ denote the symmetric Cauchy distribution with scale parameter $\gamma>0$; this distribution has density
$$
\pi(c|\gamma)=\frac1\pi\frac\gamma{\gamma^2+c^2},\qquad c\in\mathbf R.
$$
Assume a priori that $a\sim\text{Cauchy}(10)$ and $b\sim\text{Cauchy}(2.5)$ are independent (the argument for using this prior is given in a 2008 paper by Andrew Gelman and coauthors in Annals of Applied Statistics). Then construct a Metropolis-Hastings algorithm which generates a sample
from the posterior distribution of $(a,b)$.

**Hint:** You have to work with the **logarithm** of the Hastings ratio for
numerical stability. The logarithm of the prior density can be calculated as
```{r}
lprior <- function(a, b){
  dcauchy(a, scale = 10, log = TRUE) + dcauchy(b, scale = 2.5, log = TRUE)
}
```

```{r}
lposterior <- function(a, b){
  llik(a, b) + lprior(a, b)
}
myMH2 <- function(N, sigma = c(1,1), a0=1, b0=1){
  a <- b <- rep(0, N) # Empty Markov chain
  a[1] <- a0 ## initial value
  b[1] <- b0 ## initial value
  for(i in 2:N){
    a_new <- rnorm(1, a[i-1], sigma[1])
    b_new <- rnorm(1, b[i-1], sigma[2])
    logH <- lposterior(a_new, b_new) - lposterior(a[i-1], b[i-1])
    if(log(runif(1))<logH){
      a[i] <- a_new
      b[i] <- b_new
    }else{
      a[i] <- a[i-1]
      b[i] <- b[i-1]
    }
  }
  # Finally return the results
  return(mcmc(cbind(a,b)))
}
# Try it out:
chain <- myMH2(100000, sigma = c(4,4), a0 = 60, b0 = -30)
chain <- window(chain, start = 1000)
summary(chain)
plot(chain)
```

Let us plot the data (with the dose ($x$) values spread out a bit, so all the
measurements aren't on top of each other) and overlay the fitted probability of
survival for the posterior medians of the parameters (rounded to $a=58$, $b=-33$):
```{r}
plot(jitter(beet_x), beet_y, xlim = c(1.65, 1.9))
f <- function(x,a,b){ exp(a+b*x)/(1+exp(a+b*x))}
curve(f(x, a = 58, b = -33), add = TRUE, col = 2)
```

Also we can add the observed survival fractions in percent:

```{r}
tab <- table(beetles)
x <- as.numeric(rownames(tab))
y <- round(100*tab[,2]/rowSums(tab))
plot(jitter(beet_x), beet_y, xlim = c(1.65, 1.9))
curve(f(x, a = 58, b = -33), add = TRUE, col = 2)
points(x, y/100, cex = 2, pch = 3, col = 4)
text(x, y/100, labels = paste(y, "%"), pos = 2)
```


## Exercise 3: The banana

Construct a Metropolis-Hastings algorithm which has the following
(unnormalised) invariant two-dimensional density
$$
\pi(x,y) \propto \exp(-x^2/200 -(y + 0.1*x^2 -10)^2)\qquad x,y\in\mathbf R. 
$$
Notice that $\pi(x,y)$ is symmetric in $x$. Do your simulations
reproduce this symmetry?

```{r}
lf <- function(x, y){
  -x^2/200 - (y + 0.1*x^2 - 10)^2
}
myMH3 <- function(N, sigma = c(1,1), x0=0, y0=0){
  x <- y <- rep(0, N) # Empty Markov chain
  x[1] <- x0 ## initial value
  y[1] <- y0 ## initial value
  for(i in 2:N){
    x_new <- rnorm(1, x[i-1], sigma[1])
    y_new <- rnorm(1, y[i-1], sigma[2])
    logH <- lf(x_new, y_new) - lf(x[i-1], y[i-1])
    if(log(runif(1))<logH){
      x[i] <- x_new
      y[i] <- y_new
    }else{
      x[i] <- x[i-1]
      y[i] <- y[i-1]
    }
  }
  # Finally return a list of the results
  return(cbind(x,y))
}
# Try it out:
library(coda)
chain <- myMH3(10000, sigma = c(2,2), x0 = 0, y0 = 0)
plot(chain[,1], chain[,2], ylim = c(min(chain[,2]), 20), col = rgb(0,0,0,.1))
xx <- seq(-20, 20, by=.5)
yy <- seq(-10, 20, by=.5)
contour(xx, yy, outer(xx, yy, lf), levels = seq(-50, 0, by = 10), add = TRUE, col = 2)
```