####################
## Known variance ##
####################

## Reference prior
xx <- seq(115,145,len=199)
mean.post <- 130
var.post <- 25/2
plot(xx, dnorm(xx, mean.post, sd=sqrt(var.post)), lwd=3, col=pal(2)[2], type="l", xlab=expression(theta), ylab="Posterior density", ylim=c(0, 0.12))
qnorm(c(.025, .975), mean.post, sd=sqrt(var.post))

## Informative prior
n0 <- 1/4
nn <- 2+n0
mean.post <- n0/nn*120+2/nn*130
var.post <- 25/(2+n0)
lines(xx, dnorm(xx, mean.post, sd=sqrt(var.post)), lwd=3, col=pal(2)[1])
toplegend(c("Informative", "Reference"), lwd=3, col=pal(2), ncol=2)
qnorm(c(.025, .975), mean.post, sd=sqrt(var.post))

################
## Known mean ##
################

## Posterior densities
a.ref <- 1
b.ref <- 125/2
a.inf <- 3/2
b.inf <- 150/2
xx <- seq(0, 0.1, len=199)
matplot(xx, cbind(dgamma(xx, a.ref, b.ref), dgamma(xx, a.inf, b.inf)), lwd=3, type="l", lty=1, col=pal(2))
dsd <- function(x, a, b) dgamma(x^(-2), a, b)*2*x^(-3) ## Density for standard deviation
xx <- seq(0, 25, len=199)
matplot(xx, cbind(dsd(xx, a.ref, b.ref), dsd(xx, a.inf, b.inf)), lwd=3, type="l", lty=1, col=pal(2), xlab=expression(sigma), ylab="Posterior density")
toplegend(c("Reference", "Informative"), lwd=3, col=pal(2), ncol=2)

## Posterior means, modes, SDs
tau.ref <- rgamma(100000, a.ref, b.ref)
tau.inf <- rgamma(100000, a.inf, b.inf)
mean(1/sqrt(tau.ref))
mean(1/sqrt(tau.inf))
sd(1/sqrt(tau.ref))
sd(1/sqrt(tau.inf))
optimize(dsd, c(5,10), maximum=TRUE, a=a.ref, b=b.ref)$max
optimize(dsd, c(5,10), maximum=TRUE, a=a.inf, b=b.inf)$max

rev(1/sqrt(qgamma(c(.025, .975), a.ref, b.ref)))
rev(1/sqrt(qgamma(c(.025, .975), a.inf, b.inf)))
require(coda)
HPDinterval(as.mcmc(1/sqrt(tau.ref)))
HPDinterval(as.mcmc(1/sqrt(tau.inf)))

## ################
## Gibbs sampler ##
###################
source("http://web.as.uky.edu/statistics/users/pbreheny/701/S13/notes/fun.R")
gibbs <- function(x, mu, omega, n0, RSS0, B=10000) {
  n <- length(x)
  ybar <- mean(x)
  tau.init <- 1/var(x)
  theta <- matrix(NA, nrow=B, ncol=2)
  for (b in 1:B) {
    inf <- if (b==1) n*tau.init else n*theta[b-1,2]
    Mu <- (omega*mu + inf*ybar)/(omega+inf)
    Omega <- omega+inf
    theta[b,1] <- rnorm(1, Mu, sd=1/sqrt(Omega))
    RSS <- sum((x-theta[b,1])^2)
    theta[b,2] <- rgamma(1, (n0+n)/2, (RSS+RSS0)/2)
    displayProgressBar(b,B)
  }
  theta
}
x <- c(125, 130, 137)
theta <- gibbs(x, mu=0, omega=0, n0=0, RSS0=0)

## Theta
dnplot(theta[,1], xlim=c(100, 160), ylim=c(0, 0.12), xlab=expression(theta))
xx <- seq(100, 160, len=199)
lines(xx, dnorm(xx, mean=mean(x), sd=sd(x)/sqrt(length(x))), col="gray", lwd=3)
toplegend(c(expression(sigma*" unknown"), expression(sigma==6)), lwd=3, ncol=2, col=c(pal(2)[1], "gray"))
mean(theta[,1])
sd(theta[,1])
quantile(theta[,1], c(.025,.975))
mean(x)
sd(x)/sqrt(length(x))
mean(x)+c(-1,1)*qnorm(.975)*sd(x)/sqrt(length(x))

## Sigma
SD <- 1/sqrt(theta[,2])
dnplot(SD, xlim=c(0, 30), xlab=expression(sigma), ylim=c(0, 0.2))
xx <- seq(0, 30, len=199)
lines(xx, dsd(xx, 3/2, (length(x)-1)*var(x)/2), lwd=3, col="gray")
toplegend(c(expression(mu*" unknown"), expression(mu==130.6)), lwd=3, ncol=2, col=c(pal(2)[1], "gray"))
SDk <- 1/sqrt(rgamma(10000, 3/2, (length(x)-1)*var(x)/2))
mean(SD)
sd(SD)
HPDinterval(as.mcmc(SD))
mean(SDk)
sd(SDk)
HPDinterval(as.mcmc(SDk))