####Install nclbayes package:
library(nclbayes)

####Plotting partition:
par(mfrow = c(2, 2))

####Prior specs:
b = 5.41; c = 0.25; g = 2.5; h = 0.1

####Data values:
n = 23; xbar = 5.4848; s = 0.1882

####Posterior specs:
B = 5.484; C = 23.25; G = 14; H = 0.5080

####Plots for mu:
##X-axis:
x.mu = seq(4, 7, 0.001)
##Y-axis - Prior density:
prior.mu = dgt(x.mu, 2*g, b, h/(g*c))
##Y-axis - Posterior density:
posterior.mu = dgt(x.mu, 2*G, B, H/(G*C))
##Plot:
plot(x.mu, posterior.mu, type="l", lwd = 2, xlab=expression(mu), ylab="density", col="blue")
lines(x.mu, prior.mu, type="l", lty = 2, col="red")

####Plots for tau:
##X-axis:
x.tau = seq(0, 70, 0.001)
##Y-axis - Prior density:
prior.tau = dgamma(x.tau, g, h)
##Y-axis - Posterior density:
posterior.tau = dgamma(x.tau, G, H)
##Plot:
plot(x.tau, posterior.tau, type="l", lwd = 2, xlab=expression(tau), ylab="density", col="blue")
lines(x.tau, prior.tau, type="l", lty = 2, col="red")

####Plots for sigma:
##X-axis:
x.sigma = seq(0, 0.5, 0.0001)
##Y-axis - Prior density:
prior.sigma = dinvchi(x.sigma, g, h)
##Y-axis - Posterior density:
posterior.sigma = dinvchi(x.sigma, G, H)
##Plot:
plot(x.sigma, posterior.sigma, type="l", lwd = 2, xlab=expression(sigma), ylab="density", col="blue")
lines(x.sigma, prior.sigma, type="l", lty = 2, col="red")

####Probabilities:
##Actual mean density of the earth is 5.515g/cm3.  Prob that the (observed) mean density is within 0.1 of this value?
pgt(5.615, 2*G, B, H/(G*C)) - pgt(5.415, 2*G, B, H/(G*C))
##Without the data, using the prior alone:
pgt(5.615, 2*g, b, h/(g*c)) - pgt(5.415, 2*g, b, h/(g*c))


####Plots of the joint density:
##Means:
prior.mu.mean = b
posterior.mu.mean = B
prior.tau.mean = g/h
posterior.tau.mean = G/H
##Plot:
mu = seq(4.5, 6.5, len=1000)
tau = seq(0, 71, len=1000)
NGacontour(mu, tau, b, c, g, h, lty=3, col="red")
abline(h = prior.tau.mean, col="red")
abline(v = prior.mu.mean, col="red")
#NGacontour(mu, tau, b, c, g, h, lty=3, col="red", add=TRUE)

NGacontour(mu, tau, B, C, G, H, col="blue", add=TRUE)
abline(h = posterior.tau.mean, col="blue")
abline(v = posterior.mu.mean, col="blue")
#NGacontour(mu, tau, B, C, G, H, col="blue", add=TRUE)


####Confidence intervals/regions:
##Marginal prior and posterior confidence levels (Table 2.1):
##mu:
c(qgt(0.025, 2*g, b, h/(g*c)), qgt(0.975, 2*g, b, h/(g*c)))
c(qgt(0.025, 2*G, B, H/(G*C)), qgt(0.975, 2*G, B, H/(G*C)))
##tau:
hdiGamma(0.95, g, h)				#HDI
c(qgamma(0.025, g, h), qgamma(0.975, g, h))    #Equi-tailed (incorrect)
hdiGamma(0.95, G, H)				#HDI			
c(qgamma(0.025, G, H), qgamma(0.975, G, H))	#Equi-tailed (incorrect)
##sigma:
hdiInvchi(0.95, g, h)				#HDI
c(qinvchi(0.025, g, h), qinvchi(0.975, g, h))	#Equi-tailed (incorrect)			
hdiInvchi(0.95, G, H)				#HDI
c(qinvchi(0.025, G, H), qinvchi(0.975, G, H))	#Equi-tailed (incorrect)
##Confidence regions:
mu = seq(3.5, 7.5, len=1000)
tau = seq(0, 80, len=1000)
par(mfrow=c(1,1))
NGacontour(mu, tau, b, c, g, h, p=c(0.95, 0.9, 0.8), lty=3, col="red")
NGacontour(mu, tau, B, C, G, H, p=c(0.95, 0.9, 0.8), add=TRUE, col="blue")
##Zoom in:
NGacontour(mu, tau, b, c, g, h, p=c(0.95, 0.9, 0.8), lty=3, col="red", xlim=c(5.3, 5.7), ylim=c(10, 60))
NGacontour(mu, tau, B, C, G, H, p=c(0.95, 0.9, 0.8), add=TRUE, col="blue")


##Prediction:
#Prior-predictive:
c(qgt(0.025, 2*g, b, h*(c+1)/(g*c)), qgt(0.975, 2*g, b, h*(c+1)/(g*c)))
#Posterior-predictive:
c(qgt(0.025, 2*G, B, H*(C+1)/(G*C)), qgt(0.975, 2*G, B, H*(C+1)/(G*C)))