########################################
## This file is for Statistics for Data Scientists: 
##                  Monte Carlo and MCMC Simulations
## Written by James M. Flegal, February 10, 2016
########################################

set.seed(500)

########################################
## Monte Carlo Toy Example
########################################

n <- 1000
x <- rgamma(n, 3/2, scale=1)
mean(x)
y <- 1/((x+1)*log(x+3))
est <- mean(y)
est
mcse <- sd(y) / sqrt(length(y))
interval <- est + c(-1,1)*1.96*mcse
interval

## Implementing the sequential stopping rule
eps <- 0.002
len <- diff(interval)
plotting.var <- c(est, interval)
while(len > eps){
	new.x <- rgamma(n, 3/2, scale=1)
	new.y <- 1/((new.x+1)*log(new.x+3))
	y <- cbind(y, new.y)
	est <- mean(y)
	mcse <- sd(y) / sqrt(length(y))
	interval <- est + c(-1,1)*1.96*mcse
	len <- diff(interval)
	plotting.var <- rbind(plotting.var, c(est, interval))
}

## Plotting the results
temp <- seq(1000, length(y), 1000)
plot(temp, plotting.var[,1], type="l", ylim=c(min(plotting.var), 
	max(plotting.var)), main="Estimates of the Mean", xlab="Iterations", 
	ylab="Estimate")
points(temp, plotting.var[,2], type="l", col="red")
points(temp, plotting.var[,3], type="l", col="red")
legend("topright", legend=c("CI", "Estimate"), lty=c(1,1), col=c(2,1))

## Writing out a plot
pdf("GammaPlot.pdf")
plot(temp, plotting.var[,1], type="l", ylim=c(min(plotting.var), 
	max(plotting.var)), main="Estimates of the Mean", xlab="Iterations", 
	ylab="Estimate")
points(temp, plotting.var[,2], type="l", col="red")
points(temp, plotting.var[,3], type="l", col="red")
legend("topright", legend=c("CI", "Estimate"), lty=c(1,1), col=c(2,1))
dev.off()





########################################
## Introduction to MH Samplers
########################################

## Independence Metropolis sampler with Exp(theta) proposal.

ind.chain <- function(x, n, theta = 1) {
  ## if theta = 1, then this is an iid sampler
  m <- length(x)
  x <- append(x, double(n))
  for(i in (m+1):length(x)){
    x.prime <- rexp(1, rate=theta)
    u <- exp((x[(i-1)]-x.prime)*(1-theta))
    if(runif(1) < u)
      x[i] <- x.prime
    else
      x[i] <- x[(i-1)]
  }
  return(x)
}

## Random Walk Metropolis sampler with N(0,sigma) proposal.

rw.chain <- function(x, n, sigma = 1) {
  m <- length(x)
  x <- append(x, double(n))
  for(i in (m+1):length(x)){
    x.prime <- x[(i-1)] + rnorm(1, sd = sigma)
    u <- exp((x[(i-1)]-x.prime))
    u
    if(runif(1) < u && x.prime > 0)
      x[i] <- x.prime
    else
      x[i] <- x[(i-1)]
  }
  return(x)
}

## Simulations

trial0 <- ind.chain(1, 200, 1)
trial1 <- ind.chain(1, 200, 2)
trial2 <- ind.chain(1, 200, 1/2)
rw1 <- rw.chain(1, 200, .2)
rw2 <- rw.chain(1, 200, 1)
rw3 <- rw.chain(1, 200, 5)

## Plotting

par(mfrow=c(2,3))
plot.ts(trial0, ylim=c(0,6), main="IID Draws")
plot.ts(trial1, ylim=c(0,6), main="Independence with 1/2")
plot.ts(trial2, ylim=c(0,6), main="Independence with 2")
plot.ts(rw1, ylim=c(0,6), main="Random Walk with .2")
plot.ts(rw2, ylim=c(0,6), main="Random Walk with 1")
plot.ts(rw3, ylim=c(0,6), main="Random Walk with 5")
par(mfrow=c(1,1))

## Writing out a plot
pdf("MHPlot.pdf")
par(mfrow=c(2,3))
plot.ts(trial0, ylim=c(0,6), main="IID Draws")
plot.ts(trial1, ylim=c(0,6), main="Indepdence with 1/2")
plot.ts(trial2, ylim=c(0,6), main="Indepdence with 2")
plot.ts(rw1, ylim=c(0,6), main="Random Walk with .2")
plot.ts(rw2, ylim=c(0,6), main="Random Walk with 1")
plot.ts(rw3, ylim=c(0,6), main="Random Walk with 5")
par(mfrow=c(1,1))
dev.off()








set.seed(1)

########################################
## Capture-recapture Data
########################################

captured <- c(30, 22, 29, 26, 31, 32, 35)
new.captures <- c(30, 8, 17, 7, 9, 8, 5)
total.r <- sum(new.captures)

########################################
## Gibbs Sampler
########################################

gibbs.chain <- function(n, N.start = 94, alpha.start = rep(.5,7)) {
	output <- matrix(0, nrow=n, ncol=8)
	for(i in 1:n){
		neg.binom.prob <- 1 - prod(1-alpha.start)
		N.new <- rnbinom(1, 85, neg.binom.prob) + total.r
		beta1 <- captured + .5
		beta2 <- N.new - captured + .5
		alpha.new <- rbeta(7, beta1, beta2)
		output[i,] <- c(N.new, alpha.new)
		N.start <- N.new
		alpha.start <- alpha.new	
	}
	return(output)
}

########################################
## Preliminary Simulations
########################################

trial <- gibbs.chain(1000)
plot.ts(trial[,1], main = "Trace plot for N")
for(i in 1:7){
	plot.ts(trial[,(i+1)], main = paste("Trace plot for Alpha", i))
	readline("Press <return to continue")
	}

acf(trial[,1], main = "Lag Correlation plot for N")
for(i in 1:7){
	acf(trial[,(i+1)], main = paste("Lag Correlation plot for Alpha", i))
	readline("Press <return to continue")
	}

########################################
## Simulations
########################################

sim <- gibbs.chain(10000)
N <- sim[,1]
alpha1 <- sim[,2]
hist(N, freq=F, main="Estimated Marginal Posterior for N")
hist(alpha1, freq=F, main ="Estimating Marginal Posterior for Alpha 1")

library(mcmcse)

ess(N)
ess(alpha1)

estvssamp(N)
estvssamp(alpha1)

mcse(N)
mcse.q(N, .05)
mcse.q(N, .95)

mcse(alpha1)
mcse.q(alpha1, .05)
mcse.q(alpha1, .95)

current <- sim[10000,] # start from here is you need more simulations
sim <- rbind(sim, gibbs.chain(10000, N.start = current[1], alpha.start = current[2:8]))
N.big <- sim[,1]
hist(N.big, freq=F, main="Estimated Marginal Posterior for N")

ess(N)
ess(N.big)

estvssamp(N)
estvssamp(N.big)

mcse(N)
mcse(N.big)

mcse.q(N, .05)
mcse.q(N.big, .05)
mcse.q(N, .95)
mcse.q(N.big, .95)










########################################
## AR(1) Example
########################################

# The following will provide an observation from the MC 1 step ahead
ar1 <- function(m, rho, tau) {
rho*m + rnorm(1, 0, tau)
}

# Next, we will add to this program so that we can give it a Markov chain and the result will be p observations from the Markov chain.  

ar1.gen <- function(mc, p, rho, tau, q=1) {
loc <- length(mc)
junk <- double(p)
mc <- append(mc, junk)

for(i in 1:p){
j <- i+loc-1
mc[(j+1)] <- ar1(mc[j], rho, tau)
}
return(mc)
}

########################################
## AR(1) Simulations
########################################

set.seed(20)
library(mcmcse)

tau <- 1
rho <- .95
out <- 0
eps <- 0.1
start <- 1000
r <- 1000

out <- ar1.gen(out, start, rho, tau)
MCSE <- mcse(out)$se
N <- length(out)
t <- qt(.975, (floor(sqrt(N) - 1)))
muhat <- mean(out)
check <- MCSE * t

while(eps < check) {
out <- ar1.gen(out, r, rho, tau)
MCSE <- append(MCSE, mcse(out)$se)
N <- length(out)
t <- qt(.975, (floor(sqrt(N) - 1)))
muhat <- append(muhat, mean(out))
check <- MCSE[length(MCSE)] * t
}

N <- seq(start, length(out), r) 
t <- qt(.975, (floor(sqrt(N) - 1)))
half <- MCSE * t
sigmahat <- MCSE*sqrt(N)
N <- seq(start, length(out), r) / 1000

pdf("Mean.pdf")
plot(N, muhat, main="Estimates of the Mean", xlab="Iterations (in 1000's)")
points(N, muhat, type="l", col="red")
abline(h=0, lwd=3)
legend("bottomright", legend=c("Observed", "Actual"), lty=c(1,1), col=c(2,1), lwd=c(1,3))
dev.off()

plot(N, sigmahat, main="Estimates of Sigma", xlab="Iterations (in 1000's)")
points(N, sigmahat, type="l", col="red")
abline(h=20, lwd=3)
legend("bottomright", legend=c("Observed", "Actual"), lty=c(1,1), col=c(2,1), lwd=c(1,3))

plot(N, MCSE, main="Monte Carlo Standard Error", xlab="Iterations (in 1000's)")
points(N, MCSE, type="l", col="red")
abline(h=0.05, lwd=3)
legend("topright", legend=c("Observed", "Cut-off"), lty=c(1,1), col=c(2,1), lwd=c(1,3))

plot(N, half, main="Calculated Half-Widths", xlab="Iterations (in 1000's)", ylab="Half-Width", ylim=c(0, .9))
points(N, half, type="l", col="red")
abline(h=0.1, lwd=3)
legend("topright", legend=c("Observed", "Cut-off"), lty=c(1,1), col=c(2,1), lwd=c(1,3))

pdf("IntervalWidth.pdf")
plot(N, 2*half, main="Calculated Interval Widths", xlab="Iterations (in 1000's)", ylab="Width", ylim=c(0, 1.8))
points(N, 2*half, type="l", col="red")
abline(h=0.2, lwd=3)
legend("topright", legend=c("Observed", "Cut-off"), lty=c(1,1), col=c(2,1), lwd=c(1,3))
dev.off()