##Data: n=100; xbar=15; standard deviation=4.5
##Prior: mu~N(10, 1/100); tau~Ga(3, 12)
##Prior means: E[mu]=10; E[tau]=0.25

##Plots of the priors:
par(mfrow=c(1,2))
x.mu = seq(9, 11, 0.0001)
y.mu = dnorm(x.mu, mean=10, sd=1/10)
x.tau = seq(0, 1, 0.0001)
y.tau = dgamma(x.tau, 3, 12)

plot(x.mu, y.mu, type="l", main="Prior for mu", xlab="mu", ylab="Density")
plot(x.tau, y.tau, type="l", main="Prior for tau", xlab="tau", ylab="Density")

##nclbayes package:
library(nclbayes)

#Run a Gibbs sampler for this example:
posterior = gibbsNormal(N=10000, initial=c(10, 0.25), priorparam=c(10, 1/100, 3, 12), n=100, xbar=15, s=4.5)
par(mfrow=c(1,2))
plot(posterior, main="All realisations", col=c(1:length(posterior)))
plot(posterior, main="All realisations", type="l")

#Produce more informative plots:
mcmcAnalysis(posterior, rows=2)

#Process the Gibbs output - remove burn-in and thin:
posterior2 = mcmcProcess(input=posterior, burnin=10, thin=1)

#Re-produce plots:
par(mfrow=c(1,2))
plot(posterior2, main="After deleting first 10", col=c(1:length(posterior)))
plot(posterior2, main="After deleting first 10", type="l")
mcmcAnalysis(posterior2, rows=2)

#Numerical summaries:
cor(posterior2)
mcmcCi(posterior2, level=0.95)

##What about sigma?
sigma = 1/sqrt(posterior2[ ,2])
summary(sigma)
sd(sigma)
quantile(sigma, probs=c(0.025, 0.975))

##Other graphical summaries?
par(mfrow=c(1,3))
boxplot(posterior2[,1], main="Posterior for mu", ylab="mu")
boxplot(posterior2[,2], main="Posterior for tau", ylab="tau")
boxplot(sigma, main="Posterior for sigma", ylab="sigma")

par(mfrow=c(1,2))
plot(x.mu, y.mu, type="l", lty=2, xlab="mu", ylab="Density", xlim=c(9,17), ylim=c(0,5))
lines(density(posterior2[,1]), lwd=2)
plot(x.tau, y.tau, type="l", lty=2, xlab="tau", ylab="Density", xlim=c(0,0.9), ylim=c(0,55))
lines(density(posterior2[,2]), lwd=2)