# The example given as “Listing 3.1” in Reich & Ghosh
# Also discussed in class on Tue Feb 23
# Demonstrates numerical optimization and numerical integration
# based on on approximation to 2D posterior density
# Normal(mu, sigma^2) likelihood, with mu ~ Normal and
# sigma ~ Uniform priors
# Data
rm(list=ls()); y <- c(2.68, 1.18, -0.97, -0.98, -1.03); # Hyperparameters mu0 <- 0; tau=100; # Prior mean and SD for mu # Compute posterior and sketch a contour plot mu <- seq(-10, 10, .05); I <- length(mu); I; sigma <- seq(.02, 10, .02); J <- length(sigma); J; post <- matrix(nrow=I, ncol=J); for(i in 1:I){ for(j in 1:J){ post[i,j] <- dnorm(mu[i], mu0, tau, log=T) + sum(dnorm(y, mu[i], sigma[j], log=T)) }} maxie <- max(post); maxie; post <- post - maxie; rm(maxie); post <- exp(post) contours <- c(.001, .01, seq(.05, .95, .10)) contour(mu, sigma, post, levels=contours, drawlabels=F, xlab="mu", ylab="sigma", main="p(mu,sigma|y)" ) # Okay! # MAP Estimation neg.log.post <- function(theta, y) { g <- 0 - dnorm(theta[1], 0, tau, log=T) - sum(dnorm(y, theta[1], theta[2], log=T)) return(g) } inits <- c(mean(y), sd(y)) MAP <- optim(inits, neg.log.post, y=y) MAP MAP <- MAP$par points(MAP[1], MAP[2], col="blue", pch=19, cex=.75) # Compute the posterior mean post.mu <- apply(post, 1, sum) # marginal of mu post.sigma <- apply(post, 2, sum) # marginal of sigma pm.mu <- sum(mu*post.mu) / sum(post.mu); pm.mu; pm.sigma <- sum(sigma*post.sigma) / sum(post.sigma); pm.sigma; points(pm.mu, pm.sigma, col="red", cex=.75, pch=8) legend("bottomleft", inset=.05, pch=c(19,8), col=c("blue","red"), legend=c("MAP Estimate", "Posterior mean") ) # E(mu/sigma | y) h <- outer(mu, sigma, "/") sum(h * post) / sum(post)