Priors 1: Live code

Exploring informative priors

## exploring informative priors

## feel free to change the values a1 and b1 to other values of shape and rate 
a1 <- 5
b1 <- 45

ggplot()+
  stat_function(fun = dbeta, args = list(a1 , b1)) +
  xlab(expression(theta)) +
  ylab(expression(p(theta)))+
  theme(text= element_text(size = 18))

Exploring weakly informative priors

## exploring weakly informative priors

## feel free to change the values a2, b2 to other values of shape and rate 
a2 <- 2
b2 <- 2

ggplot()+
  stat_function(fun = dbeta, args = list(a2, b2)) +
  xlab(expression(theta)) +
  ylab(expression(p(theta)))+
  theme(text= element_text(size = 18))

Collect data

# observed data 
# Feel free to change these values to see the effect of strength of data!!
n <- 34
y <- 2

# visualize likelihood
theta_vec <- seq(0, 1, 0.01)
like <- dbinom(y, n,theta_vec )
ggplot(data.frame(theta = theta_vec, likelihood = dbinom(y, n,theta_vec )),
       aes(x = theta, y = likelihood)) +
  geom_line() +
  theme(text = element_text(size = 18)) +
  xlab(expression(theta)) +
  ggtitle("Likelihood of observed data")

Visualize priors

# plot priors
p1 <- ggplot() +
  stat_function(fun = dbeta, args = list(a1, b1), 
                col = "blue") +
  theme(text = element_text(size = 18)) +
  xlab(expression(theta)) +
  ylab(expression(p(theta))) +
  ggtitle("Informative prior")

p2 <- ggplot() +
  stat_function(fun = dbeta, args = list(a2, b2), 
                col = "blue") +
  theme(text = element_text(size = 18)) +
  xlab(expression(theta)) +
  ylab(expression(p(theta))) +
  ggtitle("Weakly informative prior")

ggarrange(p1, p2, ncol = 2, common.legend = T, legend = "bottom")

Visualize posteriors

# constant for scaled likelihood
c1 <- dbeta( (a1 + y - 1)/(a1 + b1 +n - 2), a1 + y, b1 + n - y) / 
  dbinom(y, n, y/n)

# plot priors and posteriors together
p1 <- ggplot() +
  stat_function(fun = dbeta, args = list(a1, b1), 
                aes(col = "prior")) +
  stat_function(fun = dbeta, args = list(a1 + y, b1 + n - y),
                aes(col = "posterior") ) +
  geom_line(data = data.frame(theta = theta_vec, y = c1 * like), 
            aes(x= theta, y = y, col = "Scaled likelihood"),
            linetype = "dashed")+
  scale_color_manual("Distribution/function", values = c("orange", "blue", "black")) +
  theme(text = element_text(size = 18)) +
  xlab(expression(theta)) +
  ylab("Density") +
  ggtitle("Prior and Posterior", subtitle = "Informative prior")


c2 <- dbeta( (a2 + y - 1)/(a2 + b2 +n - 2), a2 + y, b2 + n - y) / 
  dbinom(y, n, y/n)
p2 <- ggplot() +
  stat_function(fun = dbeta, args = list(a2, b2), 
                aes(col = "prior")) +
  stat_function(fun = dbeta, args = list(a2 + y, b2 + n - y),
                aes(col = "posterior") ) +
   geom_line(data = data.frame(theta = theta_vec, y = c2 * like), 
            aes(x= theta, y = y, col = "scaled likelihood"),
            linetype = "dashed")+
  scale_color_manual("Distribution/function", values = c("orange", "blue", "black")) +
  theme(text = element_text(size = 18)) +
  xlab(expression(theta)) +
  ylab("Density") +
  ggtitle("Prior and Posterior", subtitle = "Weakly informative prior")

ggarrange(p1, p2, ncol = 2, common.legend = T, legend = "bottom")

# posterior probability that the proportion is less than 0.1 under informative prior
pbeta(0.1, y + a1, n - y + b1)

## [1] 0.7355695

# posterior probability that the proportion is less than 0.1 under weakly informative prior
pbeta(0.1, y + a2, n - y + b2)

## [1] 0.5135534