Exercise 2: Bayesian A/B testing using MCMC Metropolis-Hastings

The MCMCmetrop1R function from The MCMCpack package implements a Metropolis-Hastings routine. Translate the model you developed in Exercise 1 so that you can fit it using MCMCmetrop1R. It should give the same, or almost the same, result!

Hint: Here is some code that fits a Normal distribution using MCMCmetrop1R. This is not the model you want to fit, but you can use this code as a starting point.

library(MCMCpack)

# Just making up some example data
x <- rnorm(30, mean = 10, sd = 5)

# Defining a function calculating a number proportional to the log probability
log_prob <- function(par, x) {
  mean <- par[1]
  sd <- par[2]
  if(sd <= 0) {
    return(-Inf)
  }
  log_prior <- dunif(mean, -100, 100, log = TRUE) + dunif(sd, 0, 100, log = TRUE)
  log_lik <- sum(dnorm(x, mean, sd, log = TRUE))
  log_prior + log_lik
}

# Fitting the model and inspecting the fit
fit <- MCMCmetrop1R(log_prob, theta.init = c(0, 1), x = x, mcmc = 1000)

## 
## 
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## The Metropolis acceptance rate was 0.57867
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

hist(fit[,1])

median(fit[,1])

## [1] 11.08

quantile(fit[,1], c(0.025, 0.975))

##   2.5%  97.5% 
##  9.021 12.923

Scroll down for a solution. But try yourself first!

Solution (but this can be done in many ways)

library(MCMCpack)

# The data
x1 <- 4
x2 <- 8
n <- 16

log_prob <- function(par, x1, x2, n) {
  p1 <- par[1]
  p2 <- par[2]
  if(any(par > 1 | par < 0)) {
    return(-Inf)
  }
  log_prior <- dunif(p1, 0, 1, log = TRUE) + dunif(p2, 0, 1, log = TRUE)
  log_lik <- dbinom(x1, n, p1, log = TRUE) + dbinom(x2, n, p2, log = TRUE)
  log_prior + log_lik
}

fit <- MCMCmetrop1R(log_prob, c(0.5, 0.5), x1 = x1, x2 = x2, n = n, burnin = 1000, mcmc = 3000)

## 
## 
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
## The Metropolis acceptance rate was 0.54475
## @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

plot(fit)

p_diff <- fit[,2] - fit[,1]
hist(p_diff)

# The probability that underlying proportion of success is larger for x2 than for x1 
mean(p_diff > 0)

## [1] 0.9173