Given x=45, n=50 x~dbinom wih theta theta~dbeta (0.85,15)

Set theta

theta <- seq(from=0,to=1,length=10000)

Likelihood

lk <- dbinom(x=45,size=50,prob=theta)
plot(theta,lk)

Prior

prior <- dbeta(theta,85,15)
plot(theta,prior)

Marginal Likelihood

ML <- mean(lk*prior)
ML
## [1] 0.1087156

Posterior

Post <- (lk*prior/ML)
plot(theta,Post)

df_post <- data.frame(theta=theta,
                      Likelihood=lk,
                      Prior=prior,
                      Posterior=Post)
head(df_post)
##        theta    Likelihood         Prior Posterior
## 1 0.00000000  0.000000e+00  0.000000e+00         0
## 2 0.00010001 2.127252e-174 3.252755e-318         0
## 3 0.00020002 7.480861e-161 6.282941e-293         0
## 4 0.00030003 6.278285e-153 3.883475e-278         0
## 5 0.00040004 2.629462e-147 1.211898e-267         0
## 6 0.00050005 6.033930e-143 1.672241e-259         0
posti <- melt(df_post, id="theta")
head(posti)
##        theta   variable         value
## 1 0.00000000 Likelihood  0.000000e+00
## 2 0.00010001 Likelihood 2.127252e-174
## 3 0.00020002 Likelihood 7.480861e-161
## 4 0.00030003 Likelihood 6.278285e-153
## 5 0.00040004 Likelihood 2.629462e-147
## 6 0.00050005 Likelihood 6.033930e-143

Compute peak per curve

peak_df <- posti %>%
  group_by(variable) %>%
  summarise(
    theta_peak = theta[which.max(value)],
    y_mid = max(value) / 2
  )
peak_df
## # A tibble: 3 × 3
##   variable   theta_peak  y_mid
##   <fct>           <dbl>  <dbl>
## 1 Likelihood      0.900 0.0925
## 2 Prior           0.857 5.67  
## 3 Posterior       0.872 7.27
ggplot(posti, aes(x=theta,y=value,group=variable,colour = variable))+geom_line(size=1.5)+theme_bw()+
  scale_x_continuous(limits = c(0.5,1))+
facet_wrap(~variable,ncol=1, axes = "all",
           scales="free_y")+
   geom_vline(data = peak_df,
             aes(xintercept = theta_peak),
             linetype = "dashed",
             colour = "black",
             inherit.aes = FALSE) +
  
geom_text(data = peak_df,
          aes(x = theta_peak,
              y = y_mid,
              label = round(theta_peak, 3)),
          angle = 90,
          vjust = -0.5,
          size = 4,
          inherit.aes = FALSE)+
  theme_bw()