n <- 400

s <- 72

alpha.prior <- 4.5

beta.prior <- 25.5

theta.grid <- seq( 0, 0.35, length = 500 )

print( alpha.likelihood <- s + 1 )

print( beta.likelihood <- n - s + 1 )

print( alpha.posterior <- alpha.prior + s )

print( beta.posterior <- beta.prior + n - s )

plot( theta.grid, dbeta( theta.grid, alpha.posterior,
  beta.posterior ), type = 'l', lwd = 2, xlab = 'theta',
  ylab = 'density', col = 'red' )

lines( theta.grid, dbeta( theta.grid, alpha.likelihood,
  beta.likelihood ), lwd = 2, col = 'black' )

lines( theta.grid, dbeta( theta.grid, alpha.prior,
  beta.prior ), lwd = 2, col = 'blue' )

text( 0.05, 6.5, 'prior', col = 'blue' )

text( 0.125, 20, 'posterior', col = 'red' )

text( 0.25, 17.5, 'likelihood', col = 'black' )

##############################################################

# 95% bayesian interval for theta

c( qbeta( 0.025, alpha.posterior, beta.posterior ),
  qbeta( 0.975, alpha.posterior, beta.posterior ) )

# [1] 0.1432590 0.2153918