bern_bayes.R

data <- list(N = 5, y = c(0, 1, 1, 0, 1))

# log probability function
# Note: the below operations are not arithmetically 
# stable. Use for demo purposed only.
lp <- function(theta, d) {
  lp <- 0
  for (i in 1:d$N) {
    lp <- lp + log(theta) * d$y[i] + 
      log(1 - theta) * (1 - d$y[i])
  }
  return(lp)
}

lp_dbinom <- function(theta, d) {
  lp <- 0
  for (i in 1:length(theta)) {
    lp[i] <- sum(dbinom(d$y,
                        size = 1,
                        prob = theta[i],
                        log = TRUE))
  }
  return(lp)
}

lp(c(0.6, 0.7), data)
lp_dbinom(c(0.6, 0.7), data)

library(ggplot2)
n <- 250
theta <- seq(0.001, 0.999, length.out = n)
likelihood <- lp(theta = theta, data)
prior <- log(dbeta(theta, 1, 1))
log_posterior <- likelihood + prior
posterior <- exp(log_posterior)
posterior <- posterior / sum(posterior)
post_draws <- sample(theta,
                     size = 1e5,
                     replace = TRUE,
                     prob = posterior)
post_density <- density(post_samples)
mle <- sum(data$y) / data$N
mle
d <- data.frame(x = post_density$x, y = post_density$y)
p <- ggplot(d, aes(x, y))
p <- p + geom_line(size = 5, alpha = 1/7) +  
  geom_vline(xintercept = mle, colour = "red") + theme_bw() + 
  xlab(expression(theta)) + ylab("")
p

# conjugate posterior
theta_conj <- dbeta(theta, 1 + 3, 1 + 5 - 3)
d <- data.frame(x = theta, y = theta_conj)
p + geom_line(data = d, aes(x, y)) + theme_bw()
p