Bayesian Billiards simultation in R based on Shiny app from Jason Bryer

bayesian-billiards.r

# Problem Definition: https://priorprobability.com/2014/04/27/bayesian-billiards/
# Referenced Shiny app: http://jason.bryer.org/posts/2016-02-21/Bayes_Billiards_Shiny.html

library(dplyr)

set.seed(424242)
true_p <- runif(1)

num_draws <- 1000
draws <- sample(c(0, 1), num_draws, replace=TRUE, prob=c(1-true_p, true_p))

# Summarizes Beta posterior with (default) Uniform prior
summarize_posterior <- function(row, prior_a=1, prior_b=1, num_samples=10000) {
  post_samples <- rbeta(num_samples,
                        prior_a + row$successes,
                        prior_b + row$failures)
  data_frame(
    successes=row$successes,
    failures=row$failures,
    min=min(post_samples),
    `2.5%`=quantile(post_samples, probs=.025),
    mean=mean(post_samples),
    median=median(post_samples),
    `97.5%`=quantile(post_samples, probs=.975),
    max=max(post_samples)
  )
}

draw_table <- data_frame(successes=cumsum(draws),
                         failures=seq_len(num_draws) - successes) %>%
  rowwise() %>%
  do(summarize_posterior(.))

# How far from true value is the posterior mean?
bias <- draw_table$mean - true_p
plot(bias, type="l")
abline(h=0, col="red")

# How wide are the intervals as we update our posterior?
interval_widths <- draw_table %>%
  summarize(interval_width=`97.5%` - `2.5%`)
plot(unlist(interval_widths), type="l")