## code to accompany the MC example from lecture 5
## of PSTAT 215A, Bayesian Inference, UCSB
##
## Robert B. Gramacy adapted from Peter Hoff

## ECDF example with S=10
## par(mfrow=c(1,2)) 
S <- 10
theta <- rgamma(S,68,45)
plot(ecdf(theta), verticals=TRUE, bty="n", lwd=2,
     ylab="F-hat(x)", xlab="theta", main="")
x <- seq(0.5,2.5, length=1000)
lines(x, pgamma(x, 68, 45), col=2, lty=2, lwd=2)

## the corresponding histogram and kernel density
hist(theta, xlim=range(x), freq=FALSE, main="",
     ylab="density", bty="n")
lines(density(theta), lwd=2, xlim=range(x))
lines(x, dgamma(x, 68, 45), col=2, lty=2, lwd=2)
legend("topright", c("kernel", "truth"), col=1:2,
       lty=1:2, lwd=2, bty="n")

## try for larger S
## S <- 100
## S <- 1000

## returning the the education/birthrates example
a <- 2; b <- 1  ## prior
sy2 <- 66; n2 <- 44 ## data

## gether some MC samples from the posterior
theta.mc10 <- rgamma(10, a+sy2, b+n2)
theta.mc100 <- rgamma(100, a+sy2, b+n2)
theta.mc1000 <- rgamma(1000, a+sy2, b+n2)
theta.mc10000 <- rgamma(10000, a+sy2, b+n2)

## mean comparison

## truth
m.truth <- (a + sy2)/(b+n2)
m.truth

## look at the MC estimates
mean(theta.mc10)
mean(theta.mc100)
mean(theta.mc1000)
mean(theta.mc10000)

## P(theta < 1.75) comparison

## truth
p.truth <- pgamma(1.75,68,45)
p.truth

## look at the MC estimates
mean(theta.mc10 < 1.75)
mean(theta.mc100 < 1.75)
mean(theta.mc1000 < 1.75)
mean(theta.mc10000 < 1.75)

## 95% CI comparison

## truth
q.truth <- qgamma(c(0.025, 0.975),68,45)

## look at the MC estimates
quantile(theta.mc10, c(0.025, 0.975))
quantile(theta.mc100, c(0.025, 0.975))
quantile(theta.mc1000, c(0.025, 0.975))
quantile(theta.mc10000, c(0.025, 0.975))


## plotting convergence diagnostics
theta <- theta.mc1000
p <- q1 <- q2 <- m <- rep(NA, length(theta))
for(i in 1:length(theta)) {
  m[i] <- mean(theta[1:i])
  p[i] <- mean(theta[1:i] < 1.75)
  q1[i] <- quantile(theta[1:i], 0.025)
  q2[i] <- quantile(theta[1:i], 0.975)
}

## convergence of mean
plot(m, type="l", ylab="cumulative mean",
     xlab="# of MC samples", lwd=2, bty="n")
abline(h=m.truth, lty=2, col=2, lwd=2)

## convergence of P(theta < 1.75)
plot(p, type="l", ylab="cumulative cdf of 1.75",
     xlab="# of MC samples", lwd=2, bty="n")
abline(h=p.truth, lty=2, col=2, lwd=2)

## convergence of quantiles
plot(q1, type="l", ylim=range(c(q1,q2)), ylab="cumulative CI",
     xlab="# of MC samples", lwd=2, bty="n")
lines(q2, lwd=2)
abline(h=q.truth[1], lty=2, col=2, lwd=2)
abline(h=q.truth[2], lty=2, col=2, lwd=2)


## log odds example

## prior
a <- b <- 1
S <- 10000
theta.prior.mc <- rbeta(S, a, b)
gamma.prior.mc <- log(theta.prior.mc/(1-theta.prior.mc))

## posterior
n1 <- 441
n0 <- 860 - n1
theta.post.mc <- rbeta(S, a+n1, b+n0)
gamma.post.mc <- log(theta.post.mc/(1-theta.post.mc))

## calculate kernel density estimates
gamma.dprior <- density(gamma.prior.mc)
gamma.dpost <- density(gamma.post.mc)

## plot the prior and posterior kernel densities
ry <- range(c(gamma.dprior$y, gamma.dpost$y))
rx <- range(c(gamma.dprior$x, gamma.dpost$x))
plot(gamma.dpost, main="", ylim=ry, ylab="density",
     xlim=rx, xlab="gamma", lwd=2, bty="n")
lines(gamma.dprior, lwd=2, col=2, lty=2)