## Set the paths for cache and figure library(methods) library(knitr) basename <- gsub(".Rmd", "", knitr:::knit_concord$get('infile')) opts_chunk$set(fig.path = paste("figure/", basename, "-", sep=""), cache.path = paste("cache/", basename, "/", sep="")) opts_chunk$set(cache = 1) opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE, comment = NA, verbose = TRUE)
library(nonparametricbayes) models <- c("Myers","Allen") parameters <- list(Myers = list( c(r=1.5 + rnorm(1, 0, .1), theta=2.5 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .2)), c(r=1.5 + rnorm(1, 0, .1), theta=2.5 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .2)), c(r=1.5 + rnorm(1, 0, .1), theta=2.5 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .2))), Allen = list( c(r=2 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .1), C=5 + rnorm(1, 0, .2)), c(r=2 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .1), C=5 + rnorm(1, 0, .2)), c(r=2 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .1), C=5 + rnorm(1, 0, .2))) ) nuisance_pars <- c("sigma_g") nuisance_values <- list(sigma_g = c(0.01, 0.05, 0.1)) replicates <- c(1111, 2222, 3333, 4444, 5555, 6666, 7777, 8888) # seeds ## @knitr posterior-mode require(modeest) posterior.mode <- function(x) { mlv(x, method="shorth")$M }
sensitivity <- function(model, parameters, nuisance, seed){ if(model == "Myers") f <- Myers else if(model == "Allen") f <- RickerAllee sigma_g <- nuisance[["sigma_g"]] z_g <- function() rlnorm(1, 0, sigma_g) p <- parameters f = f p = p z_g = z_g x_grid = seq(0, 15, length=50) h_grid = x_grid profit = function(x,h) pmin(x, h) delta = 0.01 OptTime = 50 # stationarity with unstable models is tricky thing reward = 0 xT = 0 Xo = 5.5 # observations start from x0 = 8 # simulation under policy starts from Tobs = 40 MaxT = 1000 # timeout for value iteration convergence # replicate over random seed yields <- sapply(seed, function(seed_i){ set.seed(seed_i) ## Simulate data x <- numeric(Tobs) x[1] <- Xo nz <- 1 for(t in 1:(Tobs-1)) x[t+1] = z_g() * f(x[t], h=0, p=p) X = c(rep(0,nz), pmax(rep(0,Tobs-1), x[1:(Tobs-1)])) Y = c(rep(0,nz), x[2:Tobs]) ## GP Stuff ## @knitr gp-priors s2.p <- c(5,5) d.p = c(10, 1/0.1) ## @knitr gp gp <- gp_mcmc(X, y=Y, n=1e5, s2.p = s2.p, d.p = d.p) gp_dat <- gp_predict(gp, x_grid, burnin=1e4, thin=300) matrices_gp <- gp_transition_matrix(gp_dat$Ef_posterior, gp_dat$Vf_posterior, x_grid, h_grid) opt_gp <- value_iteration(matrices_gp, x_grid, h_grid, MaxT, xT, profit, delta, reward) ## @knitr jags-setup y <- x N <- length(x); jags.data <- list("N","y") n.chains <- 6 n.iter <- 1e6 n.burnin <- floor(10000) n.thin <- max(1, floor(n.chains * (n.iter - n.burnin)/1000)) n.update <- 10 ## @knitr common-priors stdQ_prior_p <- c(1e-6, 100) stdR_prior_p <- c(1e-6, .1) stdQ_prior <- function(x) dunif(x, stdQ_prior_p[1], stdQ_prior_p[2]) stdR_prior <- function(x) dunif(x, stdR_prior_p[1], stdR_prior_p[2]) ## @knitr allen-model K_prior_p <- c(0.01, 20.0) r0_prior_p <- c(0.01, 6.0) theta_prior_p <- c(0.01, 20.0) bugs.model <- paste(sprintf( "model{ K ~ dunif(%s, %s) r0 ~ dunif(%s, %s) theta ~ dunif(%s, %s) stdQ ~ dunif(%s, %s)", K_prior_p[1], K_prior_p[2], r0_prior_p[1], r0_prior_p[2], theta_prior_p[1], theta_prior_p[2], stdQ_prior_p[1], stdQ_prior_p[2]), " iQ <- 1 / (stdQ * stdQ); y[1] ~ dunif(0, 10) for(t in 1:(N-1)){ mu[t] <- log(y[t]) + r0 * (1 - y[t]/K)* (y[t] - theta) / K y[t+1] ~ dlnorm(mu[t], iQ) } }") writeLines(bugs.model, "allen_process.bugs") ## @knitr allen-priors K_prior <- function(x) dunif(x, K_prior_p[1], K_prior_p[2]) r0_prior <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2]) theta_prior <- function(x) dunif(x, theta_prior_p[1], theta_prior_p[2]) par_priors <- list(K = K_prior, deviance = function(x) 0 * x, r0 = r0_prior, theta = theta_prior, stdQ = stdQ_prior) ## @knitr allen-mcmc jags.params=c("K","r0","theta","stdQ") # be sensible about the order here jags.inits <- function(){ list("K"= 10 * rlnorm(1,0, 0.1), "r0"= 1 * rlnorm(1,0, 0.1) , "theta"= 5 * rlnorm(1,0, 0.1) , "stdQ"= abs( 0.1 * rlnorm(1,0, 0.1)), .RNG.name="base::Wichmann-Hill", .RNG.seed=123) } set.seed(1234) # parallel refuses to take variables as arguments (e.g. n.iter = 1e5 works, but n.iter = n doesn't) allen_jags <- do.call(jags, list(data=jags.data, inits=jags.inits, jags.params, n.chains=n.chains, n.iter=n.iter, n.thin=n.thin, n.burnin=n.burnin, model.file="allen_process.bugs")) # Run again iteratively if we haven't met the Gelman-Rubin convergence criterion recompile(allen_jags) # required for parallel allen_jags <- do.call(autojags, list(object=allen_jags, n.update=n.update, n.iter=n.iter, n.thin = n.thin)) ## @knitr allen-traces tmp <- lapply(as.mcmc(allen_jags), as.matrix) # strip classes the hard way... allen_posteriors <- melt(tmp, id = colnames(tmp[[1]])) ## @knitr allen-output A <- allen_posteriors A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index pardist <- acast(A, index ~ variable) ## @knitr ricker-model K_prior_p <- c(0.01, 40.0) r0_prior_p <- c(0.01, 20.0) bugs.model <- paste(sprintf( "model{ K ~ dunif(%s, %s) r0 ~ dunif(%s, %s) stdQ ~ dunif(%s, %s)", K_prior_p[1], K_prior_p[2], r0_prior_p[1], r0_prior_p[2], stdQ_prior_p[1], stdQ_prior_p[2]), " iQ <- 1 / (stdQ * stdQ); y[1] ~ dunif(0, 10) for(t in 1:(N-1)){ mu[t] <- log(y[t]) + r0 * (1 - y[t]/K) y[t+1] ~ dlnorm(mu[t], iQ) } }") writeLines(bugs.model, "ricker_process.bugs") ## @knitr ricker-priors K_prior <- function(x) dunif(x, K_prior_p[1], K_prior_p[2]) r0_prior <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2]) par_priors <- list(K = K_prior, deviance = function(x) 0 * x, r0 = r0_prior, stdQ = stdQ_prior) ## @knitr ricker-mcmc jags.params=c("K","r0", "stdQ") jags.inits <- function(){ list("K"= 10 * rlnorm(1,0,.5), "r0"= rlnorm(1,0,.5), "stdQ"=sqrt(0.05) * rlnorm(1,0,.5), .RNG.name="base::Wichmann-Hill", .RNG.seed=123) } set.seed(12345) ricker_jags <- do.call(jags, list(data=jags.data, inits=jags.inits, jags.params, n.chains=n.chains, n.iter=n.iter, n.thin=n.thin, n.burnin=n.burnin, model.file="ricker_process.bugs")) recompile(ricker_jags) ricker_jags <- do.call(autojags, list(object=ricker_jags, n.update=n.update, n.iter=n.iter, n.thin = n.thin, progress.bar="none")) ## @knitr ricker-traces tmp <- lapply(as.mcmc(ricker_jags), as.matrix) # strip classes the hard way... ricker_posteriors <- melt(tmp, id = colnames(tmp[[1]])) names(ricker_posteriors) = c("index", "variable", "value", "chain") ## @knitr ricker-output A <- ricker_posteriors A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index ricker_pardist <- acast(A, index ~ variable) ## @knitr myers-model r0_prior_p <- c(.0001, 10.0) theta_prior_p <- c(.0001, 10.0) K_prior_p <- c(.0001, 40.0) bugs.model <- paste(sprintf( "model{ r0 ~ dunif(%s, %s) theta ~ dunif(%s, %s) K ~ dunif(%s, %s) stdQ ~ dunif(%s, %s)", r0_prior_p[1], r0_prior_p[2], theta_prior_p[1], theta_prior_p[2], K_prior_p[1], K_prior_p[2], stdQ_prior_p[1], stdQ_prior_p[2]), " iQ <- 1 / (stdQ * stdQ); y[1] ~ dunif(0, 10) for(t in 1:(N-1)){ mu[t] <- log(r0) + theta * log(y[t]) - log(1 + pow(abs(y[t]), theta) / K) y[t+1] ~ dlnorm(mu[t], iQ) } }") writeLines(bugs.model, "myers_process.bugs") ## @knitr myers-priors K_prior <- function(x) dunif(x, K_prior_p[1], K_prior_p[2]) r_prior <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2]) theta_prior <- function(x) dunif(x, theta_prior_p[1], theta_prior_p[2]) par_priors <- list( deviance = function(x) 0 * x, K = K_prior, r0 = r_prior, theta = theta_prior, stdQ = stdQ_prior) ## @knitr myers-mcmc jags.params=c("r0", "theta", "K", "stdQ") jags.inits <- function(){ list("r0"= 1 * rlnorm(1,0,.1), "K"= 10 * rlnorm(1,0,.1), "theta" = 1 * rlnorm(1,0,.1), "stdQ"= sqrt(0.2) * rlnorm(1,0,.1), .RNG.name="base::Wichmann-Hill", .RNG.seed=123) } set.seed(12345) myers_jags <- do.call(jags, list(data=jags.data, inits=jags.inits, jags.params, n.chains=n.chains, n.iter=n.iter, n.thin=n.thin, n.burnin=n.burnin, model.file="myers_process.bugs")) recompile(myers_jags) myers_jags <- do.call(autojags, list(myers_jags, n.update=n.update, n.iter=n.iter, n.thin = n.thin, progress.bar="none")) ## @knitr myers-traces tmp <- lapply(as.mcmc(myers_jags), as.matrix) # strip classes myers_posteriors <- melt(tmp, id = colnames(tmp[[1]])) names(myers_posteriors) = c("index", "variable", "value", "chain") ## @knitr myers-output A <- myers_posteriors A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index myers_pardist <- acast(A, index ~ variable) allen_deviance <- - posterior.mode(pardist[,'deviance']) ricker_deviance <- - posterior.mode(ricker_pardist[,'deviance']) myers_deviance <- - posterior.mode(myers_pardist[,'deviance']) true_deviance <- 2*estf(c(p, sigma_g)) mle_deviance <- 2*estf(c(est$p, est$sigma_g)) aictable <- data.frame(Allen = allen_deviance + 2*(1+length(bayes_pars)), # +1 for noise parameter Ricker = ricker_deviance + 2*(1+length(ricker_bayes_pars)), Myers = myers_deviance + 2*(1+length(myers_bayes_pars)), row.names = c("AIC")) bictable <- data.frame(Allen = allen_deviance + log(length(x))*(1+length(bayes_pars)), Ricker = ricker_deviance + log(length(x))*(1+length(ricker_bayes_pars)), Myers = myers_deviance + log(length(x))*(1+length(myers_bayes_pars)), row.names = c("BIC")) df <- rbind(dictable, aictable, bictable) }) dat <- melt(yields, id=names(yields[[1]])) }
model <- "Allen" allen1.01 <- sensitivity(model, parameters = parameters[[model]][[1]], nuisance = c(sigma_g = nuisance_values$sigma_g[1]), seed=c(1234, 2222, 3333))
model <- "Allen" allen2.01 <- sensitivity(model, parameters = parameters[[model]][[2]], nuisance = c(sigma_g = nuisance_values$sigma_g[1]), seed=c(1234, 2222, 3333))
model <- "Allen" allen1.05 <- sensitivity(model, parameters = parameters[[model]][[1]], nuisance = c(sigma_g = nuisance_values$sigma_g[2]), seed=c(1234, 2222, 3333))
model <- "Allen" allen2.05 <- sensitivity(model, parameters = parameters[[model]][[2]], nuisance = c(sigma_g = nuisance_values$sigma_g[2]), seed=c(1234, 2222, 3333))
model <- "Myers" Myers1.01 <- sensitivity(model, parameters = parameters[[model]][[1]], nuisance = c(sigma_g = nuisance_values$sigma_g[1]), seed=c(1234, 2222, 3333))
model <- "Myers" Myers2.01 <- sensitivity(model, parameters = parameters[[model]][[2]], nuisance = c(sigma_g = nuisance_values$sigma_g[1]), seed=c(1234, 2222, 3333))
model <- "Myers" Myers1.05 <- sensitivity(model, parameters = parameters[[model]][[1]], nuisance = c(sigma_g = nuisance_values$sigma_g[2]), seed=c(1234, 2222, 3333))
model <- "Myers" Myers2.05 <- sensitivity(model, parameters = parameters[[model]][[2]], nuisance = c(sigma_g = nuisance_values$sigma_g[2]), seed=c(1234, 2222, 3333))
## Assemble into data.frame allen_dat <- rbind(allen1.01, allen1.05, allen2.01, allen2.05) myers_dat <- rbind(Myers1.01, Myers1.05, Myers2.01, Myers2.05)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.