# Example 1: Election polling

# prior: beta(alpha, beta)

alpha.p = 225

beta.p = 185

# prior information: Ask about upper 95th percentile and lower 5th percentile.

up.p = 0.59
low.p = 0.51

x = seq(0.45,0.65,0.001)
y = dbeta(x,alpha.p,beta.p)

plot(x,y,type="l",main="prior on p",ylim=c(0,32))

prior.p.95 = pbeta(up.p,alpha.p,beta.p)
prior.p.95

prior.p.05 = pbeta(low.p,alpha.p,beta.p)
prior.p.05

# prior probability interval on p

prior.p.025 = qbeta(0.025,alpha.p,beta.p)
prior.p.025

prior.p.975 = qbeta(0.975,alpha.p,beta.p)
prior.p.975

# data: x = 621, n = 1000, assume x ~ bin(n,p)
# use Monte Carlo simulation to estimate the posterior distribution

# data

x = 621

n = 1000

# mathematical solution

post.alpha.p = alpha.p + x
post.beta.p = beta.p + n - x

post.p.mean = post.alpha.p/(post.alpha.p+post.beta.p)
post.p.mean
post.p.var = (post.alpha.p*post.beta.p)/((post.alpha.p+post.beta.p)^2*(post.alpha.p+post.beta.p+1))
post.p.var

post.p.025 = qbeta(0.025,post.alpha.p,post.beta.p)
post.p.025
post.p.975 = qbeta(0.975,post.alpha.p,post.beta.p)
post.p.975

# Monte Carlo solution

M = 10000

post.sample.p = rbeta(M,post.alpha.p,post.beta.p)

ts.plot(post.sample.p)

hist(post.sample.p, xlim=c(0.45,0.65),probability = T)

p.pm = mean(post.sample.p)
p.pm
p.pv = var(post.sample.p)
p.pv
p.pstdev = sd(post.sample.p)
p.pstdev
p.pmed = median(post.sample.p)
p.pmed
p.pCI = quantile(post.sample.p,c(0.025,0.975))
p.pCI

# plot the prior and posterior on the same plot

post.y = density(post.sample.p)

x = seq(0.45,0.65,0.001)
y = dbeta(x,alpha.p,beta.p)

plot(x,y,type="l",main="prior/posterior on p",ylim=c(0,32))
lines(post.y,type="l")