######################################################################
# STAT 7630 Bayesian Statistics 
# Peng Zeng @ Auburn University
# 09-04-2025
######################################################################

######################################################################
# IQ of a person 
######################################################################

ybar = 140
theta = seq(60, 180, length = 501)
post.norm = dnorm(ybar, theta, 15/2) * dnorm(theta, 100, 12.159)
post.norm = post.norm / (sum(post.norm) * (theta[2] - theta[1]))
post.t = dnorm(ybar, theta, 15/2) * dt((theta - 100)/6.849, 2)
post.t = post.t / (sum(post.t) * (theta[2] - theta[1]))

mean.norm = sum(post.norm * theta) * (theta[2] - theta[1])
sd.norm = sqrt(sum(post.norm * (theta - mean.norm)^2) * (theta[2] - theta[1]))
mean.t = sum(post.t * theta) * (theta[2] - theta[1])
sd.t = sqrt(sum(post.t * (theta - mean.t)^2) * (theta[2] - theta[1]))
c(mean.norm, sd.norm, mean.t, sd.t)

par(mar = c(3, 4.5, 0.5, 0.5))
plot(theta, post.norm, type = "l", col = "blue", xlab = "", ylab = "density")
lines(theta, post.t, col = "red")
lines(theta, dnorm(theta, 100, 12.159), col = "blue", lty = 2)
lines(theta, dt((theta - 100)/6.849, 2)/6.849, col = "red", lty = 2)

######################################################################
# THE END
######################################################################