# STATS 4CI3/6CI3 Winter 2019

# Examples done in class on March 18-21 2019

# Bayesian Inference for a binomial probability

# The prior density
prior <- function(theta) pi/2*sin(pi*theta)

# The data
n <- 10
x <- 2

# The posterior function (up to a normalizing constant).
posterior <- function(theta, x, n)
  100*sin(pi*theta)*theta^x*(1-theta)^(n-x)

# Generate samples using Accept-Reject

# Find the maximum of the density ratio by numeric optimization
max.ratio <- optimize(function(theta) posterior(theta, 2, 10), 
                      interval=c(0,1), maximum=TRUE)
M <- max.ratio$objective
set.seed(18032019)
N <- 1000000
Y <- runif(N) # Candidates
py <- posterior(Y, x, n)/M # the acceptance probability
U <- runif(N)
theta.obs <- Y[U<py] # the sample from the posterior

# The normalized posterior density
posterior1 <- function(th, x, n) {
  K <- integrate(function(th) posterior(th,x,n), 0, 1)$value
  posterior(th,x,n)/K
}
hist(theta.obs, breaks=100, prob=T)
lines(seq(0,1,by=0.001), posterior1(seq(0,1,by=0.001),x,n), col="red")
mean(theta.obs)
var(theta.obs)

# The true posterior mean and variance evaluated by numeric integration
integrate(function(th) th*posterior1(th, 2, 10), 0, 1)$value
integrate(function(th) th^2*posterior1(th, 2, 10), 0, 1)$value-
  integrate(function(th) th*posterior1(th, 2, 10), 0, 1)$value^2

# Using MCMC we can evaluate other characteristics more easily such as
# the posterior median and posterior quantiles for credible intervals.
median(theta.obs)
quantile(theta.obs, c(0.025,0.975))

# To find quantiles without Monte Carlo we need to solve the equation
# F(x)=p for a fixed value of p.  This is generally fairly difficult and
# requires numeric methods.  
cdf <- function(x)
  integrate(function(th) posterior(th,2,10), 0, x)$value/
      integrate(function(th) posterior(th,2,10), 0, 1)$value
optimize(function(x, p) (cdf(x)-0.5)^2, c(0,1))$minimum

c(optimize(function(x, p) (cdf(x)-0.025)^2, c(0,1))$minimum,
  optimize(function(x, p) (cdf(x)-0.975)^2, c(0,1))$minimum
)

######################################################################
# Bayesian Infernence for the normal mean and variance.
######################################################################

# We use a normal(m, v) prior for mu and an gamma(a,b) prior for 1/sigma^2
# so that the prior for sigma^2 is the inverted gamma distribution.

# We saw in class that the full conditional posterior distribution
# of mu given sigma^2 is normal and the conditional posterior for
# sigma^2 given mu is the inverted gamma. 

# We will implement the Gibbs sampler to do generate samples from
# the posterior distribution.

set.seed(180320191)
x.norm <- c(17.03, 18.45, 18.59, 18.58, 18.23, 21.78, 13.71, 17.66,
            21.93, 24.40, 14.81, 26.19, 16.36, 22.97, 26.67, 24.68)
n <- length(x.norm)
# Prior parameters
m <- 10; v <- 100
a <- 3; b <- 2

# In this example I will run the chain for 2N iterations and
# discard the first half of the chain.
N <- 10000
post.out <- matrix(NA, nrow=2*N+1, ncol=2) 

# We will take our initial values from the prior distributions
post.out[1,1] <- rnorm(1,m,v)
post.out[1,2] <- 1/rgamma(1, a, scale=b)

# Now we run the Gibbs Sampler
for (i in 1:(2*N)){
  sigma2 <- post.out[i,2]
  m.post <- (n*v*mean(x.norm)+m*sigma2)/(n*v+sigma2)
  v.post <- v*sigma2/(n*v+sigma2)
  mu <- rnorm(1, m.post, sqrt(v.post))
  a.post <- n/2+a
  b.post <- 2*b/(2+b*sum((x.norm-mu)^2))
  sigma2 <- 1/rgamma(1, a.post, scale=b.post)
  post.out[i+1,] <- c(mu, sigma2)
}
post.out <- post.out[-(1:(N+1)),]

colMeans(post.out)
library(matrixStats)
colVars(post.out)

quantile(post.out[,1], c(0.025, 0.975))
quantile(post.out[,2], c(0.025, 0.975))