#########################
# Estimating Proportion #
#########################

# Plot Beta Distribution
theta = seq(0.1, 0.99, 0.01)
par(mfrow=c(3,2))
plot(theta, dbeta(theta,0.5,0.5), type="o", ylab="f(theta)")
text(locator(1), "alpha=0.5,beta=0.5")
plot(theta, dbeta(theta,1,1), type="o", ylab="f(theta)")
text(locator(1), "alpha=1,beta=1")
plot(theta, dbeta(theta,2,4), type="o", ylab="f(theta)")
text(locator(1), "alpha=2,beta=4")
plot(theta, dbeta(theta,4,2), type="o", ylab="f(theta)")
text(locator(1), "alpha=4,beta=2")
plot(theta, dbeta(theta,4,6), type="o", ylab="f(theta)")
text(locator(1), "alpha=4,beta=6")
plot(theta, dbeta(theta,40,60), type="o", ylab="f(theta)")
text(locator(1), "alpha=40,beta=60")


# Analytical solution
theta = seq(0.1, 0.99, 0.01)
pr = dbeta(theta, 100, 100)
m = 100
n = 50
lik = theta^n*(1-theta)^(m-n)
po = pr*lik
po = po/sum(po)
par(mfrow=c(1,2))
plot(theta, pr, type="o", ylab="Prior Probability")
plot(theta, po, type="o", ylab="Posterior Probability")

# R function to obtain prior and posterior probabilities, and report posterior mean
prpo = function(alpha, beta, m, n) {
theta = seq(0.1, 0.99, 0.01)
pr = dbeta(theta, alpha, beta)
lik = theta^m*(1-theta)^(n-m)
po = pr*lik
po = po/sum(po)
dis = cbind(pr, po)
po_mean = sum(theta*po)
return(list(po_mean, dis))
}
out1 = prpo(1, 1, 8, 10)
cat("Posterior Mean is ", out1[[1]], "\n")
out2 = prpo(100, 100, 8, 10)
cat("Posterior Mean is ", out2[[1]], "\n")

theta = seq(0.1, 0.99, 0.01)
par(mfrow=c(2,2))
plot(theta, out1[[2]][,1], main="Prior", xlab="theta", type = "l")
text(locator(1), "alpha=1,beta=1")
plot(theta, out1[[2]][,2], main="Posterior", xlab="theta", type = "o")
text(locator(1), "x=8,n=10")
plot(theta, out2[[2]][,1], main="Prior", xlab="theta", type = "l")
text(locator(1), "alpha=100,beta=100")
plot(theta, out2[[2]][,2], main="Posterior", xlab="theta", type = "o")
text(locator(1), "x=8,n=10")