Monte Carlo Markov Chain algorithm (MCMC)

Metropolis Hastings algorithm

Key idea: building a reversible Markov chain with $[\theta\vert y]$ as stationnary distribution

1. Initialization $\theta^{(0)}$ an admissible initial value

2. For i in 1:nIter

   • Propose a new candidate value $\theta_c^{(i)}$ sampled from a proposal distribution $g(. \vert \theta^{(i-1)})$
   • Compute Metropolis Hastings ratio $$r_i=\class{bleu}{\frac{[y\vert \theta_c^{(i)}] [\theta_c^{(i)} ] }{[y\vert \theta^{(i-1)}] [\theta^{(i-1)} ] } } \class{rouge}{\frac{g(\theta^{(i-1)} \vert \theta_c^{(i)})}{g(\theta_c^{(i)} \vert \theta^{(i-1)})}}$$
   • Define $$\theta^{(i)} =\left \lbrace \begin{array}{l} \theta_c^{(i)} \mbox{ with probability } min( r_i,1) \cr \theta^{(i-1)} \mbox{ with probability } 1-min( r_i,1) \cr \end{array}\right.$$ `

Monte Carlo Markov Chain algorithm (MCMC)

Generic Implementation

# [Y|theta]
likelihood<- function(theta, y){
}

# [theta]
dprior <- function(theta){
}

# d(theta_c|theta_i-1)
dstep <- function(theta1, theta2){
}

# sample theta_prop from theta
rstep <- function(theta){
}

ratio_MH <- function(theta_prop, theta_c, ...){
  ratio_post <- likelihood(theta_prop) * dprior(theta_prop) /
    (likelihood(theta_c) * dprior(theta_c))
  ratio_step <- dstep(theta_prop, theta_c) / 
    dstep(theta_c, theta_prop)
  return(ratio_post * ratio_step)
}

Monte Carlo Markov Chain algorithm (MCMC)

Beta Binomial model

$$Y\vert \mathcal{B}in(N, p), \quad p\vert \mathcal{B}eta(1,1)$$

Data

$N= 10, y=3$

Propose a MCMC algorithm

Monte Carlo Markov Chain algorithm (MCMC)

Implementation for a beta Binomial model

# [Y|theta]
likelihood<- function(theta, y, N){
  if(theta$p < 0 | theta$p > 1) { 
    return(0)
  } else {
    return(dbinom(x = y, size = N, prob = theta$p))
  }
}

# [theta]
dprior <- function(theta, a.prior = 1, b.prior =1){
  dbeta(theta$p, shape1 = a.prior, shape2 = b.prior)
}

Monte Carlo Markov Chain algorithm (MCMC)

Implementation for a beta Binomial model

# d(theta_c|theta_i-1)
dstep <- function(theta1, theta2, sd.explore = 0.1){
  dnorm(theta2$p, mean = theta1$p, sd = sd.explore)
}

# sample theta_prop from theta
rstep <- function(theta, sd.explore = 0.1){
  theta_prop= list()
  theta_prop$p <- rnorm(1, mean = theta$p, sd = sd.explore)
  return(theta_prop)
}

Monte Carlo Markov Chain algorithm (MCMC)

Implementation for a beta Binomial model

# sample theta_prop from theta
ratio_MH <- function(theta_prop, theta_c, 
                     a.prior = 1, b.prior = 1, sd.explore = 0.1, 
                     y = y, N = N){
  ratio_post <- likelihood(theta_prop, y, N) * dprior(theta = theta_prop, a.prior = a.prior, b.prior = b.prior) /
    (likelihood(theta_c, y, N) * dprior(theta = theta_c, a.prior = a.prior, b.prior = b.prior))
  ratio_step <- dstep(theta1 = theta_prop, theta2 = theta_c, sd.explore = sd.explore) / 
    dstep(theta1 = theta_c, theta2 = theta_prop, sd.explore = sd.explore)
  return(ratio_post * ratio_step)
}

Monte Carlo Markov Chain algorithm (MCMC)

Implementation for a beta Binomial model

## data 
y = 3; N = 10
## init 
p_init = 0.5
## prior
a.prior = 1; b.prior = 1
## algo parameters
n_iter <- 10000; sd.explore = 0.2
##sampling
theta_sample <- tibble(p = rep(NA, n_iter))
## initialisation
theta_sample$p[1] <- p_init
## algo
for(i in 2:n_iter){
  theta_prop <- rstep(theta = theta_sample[i-1,], sd.explore = sd.explore)
  u <- ratio_MH(theta_prop = theta_prop, 
                theta_c = theta_sample[i-1,], 
                a.prior = a.prior , b.prior = b.prior, sd.explore = sd.explore, y =y, N =N)
  if( runif(1) < u ){
    theta_sample[i,] <- theta_prop
  } else {
    theta_sample[i,] <- theta_sample[i-1,]
  }
}

Monte Carlo Markov Chain algorithm (MCMC)

Implementation for a beta Binomial model

theta_sample %>%  
  rowid_to_column(var = "iter") %>%
  ggplot() + aes(x = iter, y =p) + 
  geom_line()

Monte Carlo Markov Chain algorithm (MCMC)

Implementation for a beta Binomial model

p1 <- theta_sample %>% 
  ggplot() + aes(x = p) + 
  geom_histogram(aes(y =after_stat(density)), alpha = 0.5) +
  geom_density()
p1

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Monte Carlo Markov Chain algorithm (MCMC)

Implementation for a beta Binomial model

posterior_exact <- tibble(p = seq(0,1, length.out = 10001)) %>% 
  mutate(posterior = dbeta(p, shape1 =  a.prior + y, shape2 = b.prior + N - y))
p1 + geom_line(data = posterior_exact, aes(x=p, y =posterior), col = 'red')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.