set.seed(13)
n <- 50
y.obs <- rexp(n)  ### the `real' data for this demo
print("sample moments")
print(c(mean(y.obs),var(y.obs)))

t1.obs <- mean(y.obs)
t2.obs <- sqrt(var(y.obs))
t3.obs <- median(y.obs)
t4.obs <- max(y.obs)

### choice of hyperparameters
k.0 <- v.0 <- 1  ### pretty flat!
mu.0 <- 0; sig2.0 <- 1

### update
k.n <- k.0+ n
v.n <- v.0 + n
mu.n <- (k.0*mu.0+n*mean(y.obs))/k.n
sig2.n <- (v.0*sig2.0 + (n-1)*var(y.obs) + (k.0*n/k.n)*(mean(y.obs)-mu.0)^2)/v.n
print("Posterior point estimates")
print(c(mu.n,sig2.n))

### Monte Carlo sample from posterior
NREP <- 5000
tpred <- matrix(NA, NREP, 4)
set.seed(13)
for (i in 1:NREP) {
  sig2.smp <- .5*v.n*sig2.n / rgamma(1,.5*v.n)
  th.smp <- rnorm(1, mu.n, sig2.smp/k.n)
  yrep.smp <- rnorm(n, th.smp, sig2.smp)
  tpred[i,] <- c(mean(yrep.smp), sqrt(var(yrep.smp)), 
                 median(yrep.smp), max(yrep.smp))
}     

par(mfrow=c(3,2))
hist(y.obs)

frame()

hist(tpred[,1]); abline(v=t1.obs, lty=2)
print(mean(tpred[,1]<t1.obs))

hist(tpred[,2]); abline(v=t2.obs, lty=2)
print(mean(tpred[,2]<t2.obs))

hist(tpred[,3]); abline(v=t3.obs, lty=2)
print(mean(tpred[,3]<t3.obs))

hist(tpred[,4]); abline(v=t4.obs, lty=2)
print(mean(tpred[,4]<t4.obs))