## code to accompany the Normal IQ example from 
## lecture 8 of PSTAT 215A, Bayesian Inference, UCSB
##
## Robert B. Gramacy adapted from Peter Hoff


## mse.ratio:
##
## calculate the ratio of MSEs of the MLE (e)
## and Bayes estimator (b) as a function of
## n and k0

mse.ratio <- function(n, k0)
  {
    w <- n/(k0 + n)
    mse.e <- 169/n
    mse.b <- w^2 * 169/n + (1-w)^2*144
    return(mse.b/mse.e)
  }

## plotting the MSE for various settings of k0
n <- 1:50
plot(n, mse.ratio(n, 0), ylim=c(0.45, 1.05), type="l", lwd=2,
     xlab="sample size", ylab="relative MSE", bty="n")
lines(n, mse.ratio(n, 1), col=2, lty=2, lwd=2)
lines(n, mse.ratio(n, 2), col=3, lty=3, lwd=2)
lines(n, mse.ratio(n, 3), col=4, lty=4, lwd=2)
legend("center", c("k0 = 0", "k0 = 1", "k0 = 2", "k0 = 3"),
       col=1:4, lty=1:4, lwd=2, bty="n")

## smean:
##
## calculate the sample mean of mu

smean <- function(mu0, k0, n)
  {
    w <- n/(k0 + n)
    return((w*112 + (1-w)*mu0))
  }

## svar:
##
## calculate the sampling variance of mu

svar <- function(k0, n)
  {
    w <- n/(k0 + n)
     return(w^2*169/n)
  }

## comparing the sampling distribution of the
## various estimators
mu <- seq(95,130,length=1000)
plot(mu, dnorm(mu, smean(100, 0, 10), sqrt(svar(0, 10))),
     ylim=c(0,0.125), type="l", lwd=2, xlab="mean IQ",
     ylab="density", bty="n")
lines(mu, dnorm(mu, smean(100, 1, 10), sqrt(svar(1, 10))),
      col=2, lty=2, lwd=2)
lines(mu, dnorm(mu, smean(100, 2, 10), sqrt(svar(2, 10))),
      col=3, lty=3, lwd=2)
lines(mu, dnorm(mu, smean(100, 3, 10), sqrt(svar(3, 10))),
      col=4, lty=4, lwd=2)
abline(v=112, col="gray", lwd=3)
legend("right", c("k0 = 0", "k0 = 1", "k0 = 2", "k0 = 3", "truth"),
       col=c(1:4, "gray"), lty=c(1:4,1), lwd=2, bty="n")