MAP_Bootstrap.R

# Coin Flip MAP Bootstrap
library(tidyverse)
theme_set(theme_bw(base_size=18))

# Data
outcomes = c(rep(0, 5), rep(1, 15)) # True probability = 0.75
n_trials = length(outcomes)


# Beta Prior
prior_a = 1
prior_b = 1

# Sample from posterior
posterior_a = prior_a + sum(outcomes)
posterior_b = prior_b + (n_trials - sum(outcomes))
posterior_samples = rbeta(1000, posterior_a, posterior_b)

beta_mode = function(a, b){
  # Mode of Beta distribution with parameters a and b
  return((a - 1) / (a + b - 2))
}

(posterior_mode = beta_mode(posterior_a, posterior_b))
qplot(posterior_samples, bins=10) + 
  coord_cartesian(xlim=c(0, 1)) +
  geom_vline(linetype='solid', color='red', xintercept=posterior_mode, size=2)

posterior_summary = c(posterior_mode, mean(posterior_samples), sd(posterior_samples))
round(posterior_summary, 2) # > 0.64 0.63 0.01


# Bootstrap
do_bootstrap = function(i){
  resampled_outcomes = sample(outcomes, n_trials, replace=T)
  posterior_a = prior_a + sum(resampled_outcomes)
  posterior_b = prior_b + (n_trials - sum(resampled_outcomes))
  return(beta_mode(posterior_a, posterior_b))
}
bootstrap_samples = map_dbl(1:1000, do_bootstrap)

qplot(bootstrap_samples, bins=10) + 
  coord_cartesian(xlim=c(0, 1)) +
  geom_vline(linetype='solid', color='red', xintercept=mean(bootstrap_samples), size=2)

bootstrap_summary = c(mean(bootstrap_samples), sd(bootstrap_samples))
round(bootstrap_summary, 2) # > 0.64 0.63 0.01