#####################################
# Conjugate Priors with Normal Data #
#####################################

# Example 1, Chapter 3 notes (Midge data)

# Data (in mm):
x <- c(1.64, 1.70, 1.72, 1.74, 1.82, 1.82, 1.82, 1.90, 2.08)

xbar <- mean(x); n <- length(x)

my.seq <- seq(0.5,3,length=300)


#### Case when variance is known:

# Suppose we knew sigma^2 = 0.01:

# prior parameters for normal prior on mu:

my.delta <- 1.9; my.tau <- 0.3
# "95% sure" a priori that mu is between 1.3 and 2.5, say.
# What if we adjust this?

prior.mu <- dnorm(my.seq, mean=my.delta, sd=my.tau) 
# Defining the prior

posterior.mu <- dnorm(my.seq, mean = ( (my.delta/(my.tau^2) + sum(x)/0.01)/(1/(my.tau^2) + n/0.01) ), 
                sd = sqrt( 0.01*(my.tau^2) / (0.01 + n*(my.tau^2)) )  )
# Defining the posterior

plot(my.seq, prior.mu, ylim=c(-0.7,12), yaxt='n', 
     xlab="Values of mu", ylab="Probability distribution", type='l', lty=3)
# Plotting the prior
lines(my.seq, posterior.mu)
# Plotting the posterior

### Point estimates for mu:

posterior.mean <- ( (my.delta/(my.tau^2) + sum(x)/0.01)/(1/(my.tau^2) + n/0.01) )

posterior.median <- qnorm(0.50, mean = ( (my.delta/(my.tau^2) + sum(x)/0.01)/(1/(my.tau^2) + n/0.01) ), 
                sd = sqrt( 0.01*(my.tau^2) / (0.01 + n*(my.tau^2)) ) )

posterior.mode <- my.seq[posterior.mu==max(posterior.mu)]

print(paste("posterior.mean=", round(posterior.mean,3), 
      "posterior.median=", round(posterior.median,3), "posterior.mode=", round(posterior.mode,3)))

### Interval estimate for mu:

library(TeachingDemos) # loading package in order to use the hpd function

hpd.95 <- hpd(qnorm, mean = ( (my.delta/(my.tau^2) + sum(x)/0.01)/(1/(my.tau^2) + n/0.01) ), 
                sd = sqrt( 0.01*(my.tau^2) / (0.01 + n*(my.tau^2)) ), conf=0.95)
round(hpd.95, 3)

segments(hpd.95[1], 0, hpd.95[2], 0, lwd=4)
# Showing the 95% HPD interval on the graph


#### Case when mean and variance are both unknown:

# prior parameters for inverse gamma prior on sigma^2:

my.alpha <- 18; my.beta <- .34  # Based on guess of E[sig^2] = .02, var[sig^2] = .005^2

# prior parameters for normal prior on mu:
my.delta <- 1.9; s0 <- 1  # low value of s0 indicates lack of prior knowledge

library(pscl)  # loading pscl package, to use inverse gamma distribution

# May need to install the pscl package first?
# If so, type at the command line:  install.packages("pscl", dependencies=F)
# while plugged in to the internet.

my.seq.mu <- seq(1,2.5,length=300)
my.seq.sigma.sq <- seq(0.01,0.03, length=300)

posterior.1 <- rep(0,times=300)  # setting up a dummy vector to be filled in
posterior.2 <- rep(0,times=300)  # setting up a dummy vector to be filled in
joint.posterior <- matrix(0,nrow=300,ncol=300)  # setting up a dummy matrix to be filled in

for (i in 1:300)
{
for (j in 1:300)
{
posterior.1[i] <- densigamma(my.seq.sigma.sq[i], my.alpha + n/2 - 0.5, my.beta + 0.5*( sum(x^2) - n*(xbar^2) ) )
# Defining the posterior for sigma^2

posterior.2[j] <- dnorm(my.seq.mu[j], mean=((sum(x)+my.delta*s0)/(n+s0)), sd=sqrt(my.seq.sigma.sq[i]/(n+s0)) ) 
# Defining the posterior for mu given sigma^2

joint.posterior[j,i] <- posterior.1[i] * posterior.2[j]
}
}

# Plotting the joint posterior density for mu and sigma^2:
contour(x=my.seq.mu, y=my.seq.sigma.sq, z=joint.posterior)


### Point estimates for sigma^2:

p.mean.sig.sq <- (my.beta + 0.5*(sum(x^2) - n*(xbar^2)) ) / (my.alpha + n/2 - 0.5 - 1)

p.median.sig.sq <- qigamma(0.50, my.alpha + n/2 - 0.5, my.beta + 0.5*( sum(x^2) - n*(xbar^2) ) )

print(paste("posterior.mean for sigma^2=", round(p.mean.sig.sq,4), 
      "posterior.median for sigma^2=", round(p.median.sig.sq,4) ))

#### NOT QUITE THE BEST APPROACH FOR INFERENCE ABOUT mu:

### Point estimates for mu, given the pt. estimate for sigma^2:

p.mean.mu <- ((sum(x)+my.delta*s0)/(n+s0))

p.median.mu <- qnorm(0.50, mean=((sum(x)+my.delta*s0)/(n+s0)), sd=sqrt(p.median.sig.sq/(n+s0)) )

print(paste("posterior.mean for mu=", round(p.mean.mu,4), 
      "posterior.median for mu=", round(p.median.mu,4) ))

### Marginal Interval estimate for sigma^2:


hpd.95.sig.sq <- hpd(qigamma, alpha=my.alpha + n/2 - 0.5, beta=my.beta + 0.5*( sum(x^2) - n*(xbar^2) ) )

round(hpd.95.sig.sq, 4)

### Marginal Interval estimate for mu, given the pt. estimate for sigma^2:

hpd.95.mu <- hpd(qnorm, mean=((sum(x)+my.delta*s0)/(n+s0)), sd=sqrt(p.median.sig.sq/(n+s0)) )

round(hpd.95.mu, 4)
# This underestimates the uncertainty, since it does not account for the fact that 
# the point estimate of sigma^2 is not the true sigma^2.

#### A BETTER APPROACH FOR INFERENCE ABOUT mu:

# Randomly sample many values for the posterior of sigma^2:

sig.sq.values <- rigamma(n=1000000,alpha=my.alpha + n/2 - 0.5, beta=my.beta + 0.5*( sum(x^2) - n*(xbar^2) ) )

# Randomly sample many values from the posterior of mu, GIVEN the sampled values of sigma^2 above:

mu.values <- rnorm(n=1000000,mean=((sum(x)+my.delta*s0)/(n+s0)), sd=sqrt(sig.sq.values/(n+s0)) )

# Point estimates for sigma^2 and mu:

round(median(sig.sq.values),4)

round(median(mu.values),4)

# 95% HPD interval estimates for sigma^2 and mu:

round(emp.hpd(sig.sq.values),4)

round(emp.hpd(mu.values),4)


##################################
# Vague Priors with Normal Data  #
##################################

# Again using the midge data, we let 
# p(mu) = 1 and p(sigma) = 1/sigma (noninformative priors)

# Data (in mm):
x <- c(1.64, 1.70, 1.72, 1.74, 1.82, 1.82, 1.82, 1.90, 2.08)

xbar <- mean(x); n <- length(x); my.s <- sd(x)


### An interval estimate for t:

hpd.95.t <- hpd(qt, df=n-1, conf=0.95)

### Back-transforming for mu:

hpd.95.mu.vague <- (hpd.95.t)*(my.s/sqrt(n)) + xbar
round(hpd.95.mu.vague, 3)

### An interval estimate for sigma^2:

hpd.95.sig.sq.vague <- hpd(qigamma, alpha= 0.5*(n-1), beta=0.5*(n-1)*(my.s^2) )
round(hpd.95.sig.sq.vague, 3)

# How do these intervals compare to the intervals from the conjugate analysis?