albert — Sep 25, 2013, 3:52 PM
# MATH 6480 - Fall 2013
# Learning about a Normal Distribution - part 2
# CORRECTED September 25, 2013
# Suppose sleeping times of college students are N(mu, sigma).
# Prior: mu is N(8, 0.61), sigma^2 is IG(2.8, 13.5)
mu0 <- 8; tau <- 0.61; a <- 2.8; b <- 13.5
prior <- c(mu0, tau, a, b)
library(LearnBayes)
y <- with(studentdata, WakeUp - ToSleep)
y <- y[complete.cases(y)]
# Illustrate "brute force" way of simulating from the posterior
# of (mu, sigma2)
# First construct a function that computes the logarithm of the
# posterior
# (A similar function for a noninformative prior is in normchi2post)
mynormpost <- function (theta, data.prior) {
mu <- theta[1]; sig2 <- theta[2]
data <- data.prior$data
prior <- data.prior$prior
mu0 <- prior[1]; tau <- prior[2]; a <- prior[3]; b <- prior[4]
logf <- function(y, mu, sig2) -(y - mu)^2/2/sig2 - log(sig2)/2
z <- sum(logf(data, mu, sig2))
z = z - (a + 1) * log(sig2) - b / sig2 +
dnorm(mu, mean=mu0, sd=tau, log=TRUE)
return(z)
}
# Use mycontour to compute a contour graph of posterior.
mycontour(mynormpost, c(7.1, 7.7, 1.8, 3),
list(data=y, prior=prior),
xlab="Mean", ylab="Variance")
# Use simcontour to simulate from grid.
post.sim <- simcontour(mynormpost,
c(7.1, 7.7, 1.8, 3),
list(data=y, prior=prior),
1000)
Warning: Walker's alias method used: results are different from R < 2.2.0
points(post.sim)
# Say we are interested in learning about coefficient of variation
# cv = sigma / mu
post.cv <- sqrt(post.sim$y) / post.sim$x
# Draw density plot of simulated draws
plot(density(post.cv))
