# 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)

tsplot(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 <- stdev(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 <- 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")