######################################## ## 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