plots for discussing the beta-binomial model

beta-binom.R

plot_beta = function(a, b) {
  
  par(mfrow = c(1, 1))
  curve(dbeta(x, shape1 = a, shape2 = b), col = "darkorange",
        main = "beta distribution", xlab = "x", ylab = "density", lwd = 2)
  grid()

}

plot_prior_like_post = function(a, b, x, y, plot_est = FALSE) {

  # model setups
  theta = seq(0, 1, by = 0.001)
  likelihood = (theta ^ x) * (1 - theta) ^ y
  a_post = a + x
  b_post = b + y

  # plotting setup
  par(mfrow = c(3, 1))

  # prior distribution plot
  curve(dbeta(x, shape1 = a, shape2 = b), col = "darkorange", main = "prior",
        xlab = "theta", ylab = "g(theta)", lwd = 2)
  grid()
  if (plot_est) abline(v = (a - 1) / (a + b - 2), col = "darkorange")

  # likelihood plot
  plot(theta, likelihood, type = "l", main = "likelihood",
       xlab = "theta", ylab = "L(data | theta)", lwd = 2)
  grid()
  if (plot_est) abline(v = x / (x + y))

  # posterior distribution plot
  curve(dbeta(x, shape1 = a_post, shape2 = b_post), col = "dodgerblue", main = "posterior",
        xlab = "theta", ylab = "g*(theta)", lwd = 2)
  grid()
  if (plot_est) abline(v = x / (x + y))
  if (plot_est) abline(v = (a_post - 1) / (a_post + b_post - 2), col = "dodgerblue")
  if (plot_est) abline(v = (a - 1) / (a + b - 2), col = "darkorange")

}

plot_beta(a = 5, b = 5)
plot_prior_like_post(a = 2, b = 2, x = 10, y = 40)

library(manipulate)

manipulate(plot_beta(a, b),
           a = slider(min = 0, max = 10, initial = 5, ticks = TRUE, step = 0.1),
           b = slider(min = 0, max = 10, initial = 5, ticks = TRUE, step = 0.1))

manipulate(plot_prior_like_post(a, b, x, y, plot_est),
           a = slider(min = 0, max = 100, initial = 5, ticks = TRUE, step = 5),
           b = slider(min = 0, max = 100, initial = 5, ticks = TRUE, step = 5),
           x = slider(min = 0, max = 100, initial = 5, ticks = TRUE, step = 5),
           y = slider(min = 0, max = 100, initial = 5, ticks = TRUE, step = 5),
           plot_est = checkbox(initial = FALSE))