inst/examples/aic_replicates.md
In cboettig/nonparametric-bayes: Nonparametric Bayes inference for ecological models
## 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))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Allen"
allen2.01 <- sensitivity(model,
parameters = parameters[[model]][[2]],
nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Allen"
allen1.05 <- sensitivity(model,
parameters = parameters[[model]][[1]],
nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Allen"
allen2.05 <- sensitivity(model,
parameters = parameters[[model]][[2]],
nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
seed=c(1234, 2222, 3333))
Error: system is computationally singular: reciprocal condition number =
1.82252e-16
Timing stopped at: 0.144 0.004 0.148
model <- "Myers"
Myers1.01 <- sensitivity(model,
parameters = parameters[[model]][[1]],
nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Myers"
Myers2.01 <- sensitivity(model,
parameters = parameters[[model]][[2]],
nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Myers"
Myers1.05 <- sensitivity(model,
parameters = parameters[[model]][[1]],
nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Myers"
Myers2.05 <- sensitivity(model,
parameters = parameters[[model]][[2]],
nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
## Assemble into data.frame
allen_dat <- rbind(allen1.01, allen1.05,
allen2.01, allen2.05)
Error: object 'allen1.01' not found
myers_dat <- rbind(Myers1.01, Myers1.05,
Myers2.01, Myers2.05)
Error: object 'Myers1.01' not found
cboettig/nonparametric-bayes documentation built on May 13, 2019, 2:09 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## 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))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Allen"
allen2.01 <- sensitivity(model,
parameters = parameters[[model]][[2]],
nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Allen"
allen1.05 <- sensitivity(model,
parameters = parameters[[model]][[1]],
nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Allen"
allen2.05 <- sensitivity(model,
parameters = parameters[[model]][[2]],
nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
seed=c(1234, 2222, 3333))
Error: system is computationally singular: reciprocal condition number =
1.82252e-16
Timing stopped at: 0.144 0.004 0.148
model <- "Myers"
Myers1.01 <- sensitivity(model,
parameters = parameters[[model]][[1]],
nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Myers"
Myers2.01 <- sensitivity(model,
parameters = parameters[[model]][[2]],
nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Myers"
Myers1.05 <- sensitivity(model,
parameters = parameters[[model]][[1]],
nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
model <- "Myers"
Myers2.05 <- sensitivity(model,
parameters = parameters[[model]][[2]],
nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
seed=c(1234, 2222, 3333))
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 328
Initializing model
Error: replacement has 0 rows, data has 5070
## Assemble into data.frame
allen_dat <- rbind(allen1.01, allen1.05,
allen2.01, allen2.05)
Error: object 'allen1.01' not found
myers_dat <- rbind(Myers1.01, Myers1.05,
Myers2.01, Myers2.05)
Error: object 'Myers1.01' not found
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.