# Comparison of Nonparametric Bayesian Gaussian Process estimates to standard the Parametric Bayesian approach

Plotting and knitr options, (can generally be ignored)

require(modeest)
posterior.mode <- function(x) {
mlv(x, method="shorth")$M } ### Model and parameters Uses the model derived in citet("10.1080/10236190412331335373"), of a Ricker-like growth curve with an allee effect, defined in the pdgControl package, f <- RickerAllee p <- c(2, 8, 5) K <- 10 # approx, a li'l' less allee <- 5 # approx, a li'l' less Various parameters defining noise dynamics, grid, and policy costs. sigma_g <- 0.05 sigma_m <- 0.0 z_g <- function() rlnorm(1, 0, sigma_g) z_m <- function() 1+(2*runif(1, 0, 1)-1) * sigma_m x_grid <- seq(0, 1.5 * K, 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 <- allee+.5# observations start from x0 <- K # simulation under policy starts from Tobs <- 40 ### Sample Data  set.seed(1234) #harvest <- sort(rep(seq(0, .5, length=7), 5)) 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) obs <- data.frame(x = c(rep(0,nz), pmax(rep(0,Tobs-1), x[1:(Tobs-1)])), y = c(rep(0,nz), x[2:Tobs])) raw_plot <- ggplot(data.frame(time = 1:Tobs, x=x), aes(time,x)) + geom_line() raw_plot ## Maximum Likelihood set.seed(12345) estf <- function(p){ mu <- f(obs$x,0,p)
-sum(dlnorm(obs$y, log(mu), p[4]), log=TRUE) } par <- c(p[1]*rlnorm(1,0,.4), p[2]*rlnorm(1,0,.3), p[3]*rlnorm(1,0, .3), sigma_g * rlnorm(1,0,.3)) o <- optim(par, estf, method="L", lower=c(1e-5,1e-5,1e-5,1e-5)) f_alt <- f p_alt <- c(as.numeric(o$par[1]), as.numeric(o$par[2]), as.numeric(o$par[3]))
sigma_g_alt <- as.numeric(o$par[4]) est <- list(f = f_alt, p = p_alt, sigma_g = sigma_g_alt, mloglik=o$value)

Mean predictions

true_means <- sapply(x_grid, f, 0, p)
est_means <- sapply(x_grid, est$f, 0, est$p)

## Non-parametric Bayes

#inv gamma has mean b / (a - 1) (assuming a>1) and variance b ^ 2 / ((a - 2) * (a - 1) ^ 2) (assuming a>2)
s2.p <- c(5,5)
d.p = c(10, 1/0.1)

Estimate the Gaussian Process (nonparametric Bayesian fit)

gp <- gp_mcmc(obs$x, y=obs$y, n=1e5, s2.p = s2.p, d.p = d.p)
gp_dat <- gp_predict(gp, x_grid, burnin=1e4, thin=300)

Show traces and posteriors against priors

plots <- summary_gp_mcmc(gp, burnin=1e4, thin=300)

# Summarize the GP model
tgp_dat <-
data.frame(  x = x_grid,
y = gp_dat$E_Ef, ymin = gp_dat$E_Ef - 2 * sqrt(gp_dat$E_Vf), ymax = gp_dat$E_Ef + 2 * sqrt(gp_dat$E_Vf) ) ## Parametric Bayesian Models We use the JAGS Gibbs sampler, a recent open source BUGS implementation with an R interface that works on most platforms. We initialize the usual MCMC parameters; see ?jags for details. All parametric Bayesian estimates use the following basic parameters for the JAGS MCMC: 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 We will use the same priors for process and observation noise in each model, 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]) ### Parametric Bayes of correct (Allen) model We initiate the MCMC chain (init_p) using the true values of the parameters p from the simulation. While impossible in real data, this gives the parametric Bayesian approach the best chance at succeeding. y is the timeseries (recall obs has the $$x_t$$, $$x_{t+1}$$ pairs) The actual model is defined in a model.file that contains an R function that is automatically translated into BUGS code by R2WinBUGS. The file defines the priors and the model. We write the file from R as follows: K_prior_p <- c(0.01, 40.0) logr0_prior_p <- c(-6.0, 6.0) logtheta_prior_p <- c(-6.0, 6.0) bugs.model <- paste(sprintf( "model{ K ~ dunif(%s, %s) logr0 ~ dunif(%s, %s) logtheta ~ dunif(%s, %s) stdQ ~ dunif(%s, %s)", K_prior_p[1], K_prior_p[2], logr0_prior_p[1], logr0_prior_p[2], logtheta_prior_p[1], logtheta_prior_p[2], stdQ_prior_p[1], stdQ_prior_p[2]), " iQ <- 1 / (stdQ * stdQ); r0 <- exp(logr0) theta <- exp(logtheta) y[1] ~ dunif(0, 10) for(t in 1:(N-1)){ mu[t] <- y[t] * exp(r0 * (1 - y[t]/K)* (y[t] - theta) / K ) y[t+1] ~ dnorm(mu[t], iQ) } }") writeLines(bugs.model, "allen_process.bugs") Write the priors into a list for later reference K_prior <- function(x) dunif(x, K_prior_p[1], K_prior_p[2]) logr0_prior <- function(x) dunif(x, logr0_prior_p[1], logr0_prior_p[2]) logtheta_prior <- function(x) dunif(x, logtheta_prior_p[1], logtheta_prior_p[2]) par_priors <- list(K = K_prior, deviance = function(x) 0 * x, logr0 = logr0_prior, logtheta = logtheta_prior, stdQ = stdQ_prior) We define which parameters to keep track of, and set the initial values of parameters in the transformed space used by the MCMC. We use logarithms to maintain strictly positive values of parameters where appropriate. jags.params=c("K","logr0","logtheta","stdQ") # be sensible about the order here jags.inits <- function(){ list("K"= 8 * rlnorm(1,0, 0.1), "logr0"=log(2 * rlnorm(1,0, 0.1) ), "logtheta"=log( 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.parallel, 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)) #### Convergence diagnostics for Allen model R notes: this strips classes from the mcmc.list object (so that we have list of matrices; objects that reshape2::melt can handle intelligently), and then combines chains into one array. In this array each parameter is given its value at each sample from the posterior (index) for each chain. tmp <- lapply(as.mcmc(allen_jags), as.matrix) # strip classes the hard way... allen_posteriors <- melt(tmp, id = colnames(tmp[[1]])) names(allen_posteriors) = c("index", "variable", "value", "chain") ggplot(allen_posteriors) + geom_line(aes(index, value)) + facet_wrap(~ variable, scale="free", ncol=1) allen_priors <- ddply(allen_posteriors, "variable", function(dd){ grid <- seq(min(dd$value), max(dd$value), length = 100) data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})

ggplot(allen_posteriors, aes(value)) +
stat_density(geom="path", position="identity", alpha=0.7) +
geom_line(data=allen_priors, aes(x=value, y=density), col="red") +
facet_wrap(~ variable, scale="free", ncol=3)

Reshape the posterior parameter distribution data, transform back into original space, and calculate the mean parameters and mean function

A <- allen_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
pardist <- acast(A, index ~ variable)
# pardist <- acast(allen_posteriors[2:3], 1:table(allen_posteriors$variable)[1] ~ variable) # NOT SURE WHY THIS FAILS # transform model parameters back first pardist[,"logr0"] = exp(pardist[,"logr0"]) pardist[,"logtheta"] = exp(pardist[,"logtheta"]) colnames(pardist)[colnames(pardist)=="logtheta"] = "theta" colnames(pardist)[colnames(pardist)=="logr0"] = "r0" bayes_coef <- apply(pardist,2, posterior.mode) bayes_pars <- unname(c(bayes_coef["r0"], bayes_coef["K"], bayes_coef["theta"])) # parameters formatted for f allen_f <- function(x,h,p) unname(RickerAllee(x,h, unname(p[c("r0", "K", "theta")]))) allen_means <- sapply(x_grid, f, 0, bayes_pars) bayes_pars [1] 0.01929 7.72010 0.06654 head(pardist)  K deviance r0 theta stdQ 170 21.327 45.03 0.010484 3.826031 0.3737 171 14.261 45.02 0.029587 0.019707 0.4387 172 7.828 40.50 0.129248 2.769482 0.3795 173 29.450 45.66 0.044378 0.006615 0.4461 174 37.580 44.89 0.006334 1.252886 0.3766 175 20.168 45.05 0.033854 0.144171 0.4294 ## Parametric Bayes based on the structurally wrong model (Ricker) K_prior_p <- c(0.01, 40.0) logr0_prior_p <- c(-6.0, 6.0) bugs.model <- paste(sprintf( "model{ K ~ dunif(%s, %s) logr0 ~ dunif(%s, %s) stdQ ~ dunif(%s, %s)", K_prior_p[1], K_prior_p[2], logr0_prior_p[1], logr0_prior_p[2], stdQ_prior_p[1], stdQ_prior_p[2]), " iQ <- 1 / (stdQ * stdQ); r0 <- exp(logr0) y[1] ~ dunif(0, 10) for(t in 1:(N-1)){ mu[t] <- y[t] * exp(r0 * (1 - y[t]/K) ) y[t+1] ~ dnorm(mu[t], iQ) } }") writeLines(bugs.model, "ricker_process.bugs") Compute prior curves K_prior <- function(x) dunif(x, K_prior_p[1], K_prior_p[2]) logr0_prior <- function(x) dunif(x, logr0_prior_p[1], logr0_prior_p[2]) par_priors <- list(K = K_prior, deviance = function(x) 0 * x, logr0 = logr0_prior, stdQ = stdQ_prior) We define which parameters to keep track of, and set the initial values of parameters in the transformed space used by the MCMC. We use logarithms to maintain strictly positive values of parameters where appropriate. # Uniform priors on standard deviation terms jags.params=c("K","logr0", "stdQ") jags.inits <- function(){ list("K"=10 * rlnorm(1,0,.5), "logr0"=log(1) * 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.parallel, 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) Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 251 Initializing model Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 251 Initializing model Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 251 Initializing model Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 251 Initializing model Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 251 Initializing model Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 251 Initializing model ricker_jags <- do.call(autojags, list(object=ricker_jags, n.update=n.update, n.iter=n.iter, n.thin = n.thin, progress.bar="none")) #### Convergence diagnostics for parametric bayes Ricker model 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") ggplot(ricker_posteriors) + geom_line(aes(index, value)) + facet_wrap(~ variable, scale="free", ncol=1) ricker_priors <- ddply(ricker_posteriors, "variable", function(dd){ grid <- seq(min(dd$value), max(dd$value), length = 100) data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})
# plot posterior distributions
ggplot(ricker_posteriors, aes(value)) +
stat_density(geom="path", position="identity", alpha=0.7) +
geom_line(data=ricker_priors, aes(x=value, y=density), col="red") +
facet_wrap(~ variable, scale="free", ncol=2)

Reshape posteriors data, transform back, calculate mode and corresponding function.

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)
ricker_pardist[,"logr0"] = exp(ricker_pardist[,"logr0"]) # transform model parameters back first
colnames(ricker_pardist)[colnames(ricker_pardist)=="logr0"] = "r0"
bayes_coef <- apply(ricker_pardist,2, posterior.mode) # much better estimates from mode then mean
ricker_bayes_pars <- unname(c(bayes_coef["r0"], bayes_coef["K"]))

ricker_f <- function(x,h,p){
sapply(x, function(x){
x <- pmax(0, x-h)
pmax(0, x * exp(p["r0"] * (1 - x / p["K"] )) )
})
}
ricker_means <- sapply(x_grid, Ricker, 0, ricker_bayes_pars[c(1,2)])

head(ricker_pardist)
         K deviance       r0   stdQ
170 30.871    44.59 0.011138 0.4164
171 15.920    44.30 0.003492 0.4062
172  8.630    46.04 0.040291 0.4962
173  7.935    37.59 0.172355 0.3977
174  8.622    42.76 0.096497 0.3325
175  8.279    40.48 0.143866 0.4316
ricker_bayes_pars
[1] 0.008073 8.011228

## Myers Parametric Bayes

logr0_prior_p <- c(-6.0, 6.0)
logtheta_prior_p <- c(-6.0, 6.0)
logK_prior_p <- c(-6.0, 6.0)

bugs.model <-
paste(sprintf(
"model{
logr0    ~ dunif(%s, %s)
logtheta    ~ dunif(%s, %s)
logK    ~ dunif(%s, %s)
stdQ ~ dunif(%s, %s)",
logr0_prior_p[1], logr0_prior_p[2],
logtheta_prior_p[1], logtheta_prior_p[2],
logK_prior_p[1], logK_prior_p[2],
stdQ_prior_p[1], stdQ_prior_p[2]),

"
iQ <- 1 / (stdQ * stdQ);
r0 <- exp(logr0)
theta <- exp(logtheta)
K <- exp(logK)

y[1] ~ dunif(0, 10)
for(t in 1:(N-1)){
mu[t] <- r0 * pow(abs(y[t]), theta) / (1 + pow(abs(y[t]), theta) / K)
y[t+1] ~ dnorm(mu[t], iQ)
}
}")
writeLines(bugs.model, "myers_process.bugs")
logK_prior     <- function(x) dunif(x, logK_prior_p[1], logK_prior_p[2])
logr_prior     <- function(x) dunif(x, logr0_prior_p[1], logr0_prior_p[2])
logtheta_prior <- function(x) dunif(x, logtheta_prior_p[1], logtheta_prior_p[2])
par_priors <- list( deviance = function(x) 0 * x, logK = logK_prior,
logr0 = logr_prior, logtheta = logtheta_prior,
stdQ = stdQ_prior)
jags.params=c("logr0", "logtheta", "logK", "stdQ")
jags.inits <- function(){
list("logr0"=log(rlnorm(1,0,.1)),
"logK"=log(10 * rlnorm(1,0,.1)),
"logtheta" = log(2 * rlnorm(1,0,.1)),
"stdQ"=sqrt(0.5) * rlnorm(1,0,.1),
.RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
set.seed(12345)
myers_jags <- do.call(jags.parallel,
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)
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 291

Initializing model

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 291

Initializing model

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 291

Initializing model

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 291

Initializing model

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 291

Initializing model

Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 291

Initializing model
myers_jags <- do.call(autojags,
list(myers_jags, n.update=n.update, n.iter=n.iter,
n.thin = n.thin, progress.bar="none"))

Convergence diagnostics for parametric bayes

tmp <- lapply(as.mcmc(myers_jags), as.matrix) # strip classes the hard way...
myers_posteriors <- melt(tmp, id = colnames(tmp[[1]]))
names(myers_posteriors) = c("index", "variable", "value", "chain")

ggplot(myers_posteriors) + geom_line(aes(index, value)) +
facet_wrap(~ variable, scale="free", ncol=1)
par_prior_curves <- ddply(myers_posteriors, "variable", function(dd){
grid <- seq(min(dd$value), max(dd$value), length = 100)
data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid)) }) ggplot(myers_posteriors, aes(value)) + stat_density(geom="path", position="identity", alpha=0.7) + geom_line(data=par_prior_curves, aes(x=value, y=density), col="red") + facet_wrap(~ variable, scale="free", ncol=3) 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) ### myers_pardist <- acast(myers_posteriors[2:3], 1:table(myers_posteriors$variable) ~ variable)
myers_pardist[,"logK"] = exp(myers_pardist[,"logK"]) # transform model parameters back first
myers_pardist[,"logr0"] = exp(myers_pardist[,"logr0"]) # transform model parameters back first
myers_pardist[,"logtheta"] = exp(myers_pardist[,"logtheta"]) # transform model parameters back first
colnames(myers_pardist)[colnames(myers_pardist)=="logK"] = "K"
colnames(myers_pardist)[colnames(myers_pardist)=="logr0"] = "r0"
colnames(myers_pardist)[colnames(myers_pardist)=="logtheta"] = "theta"

bayes_coef <- apply(myers_pardist,2, posterior.mode) # much better estimates
myers_bayes_pars <- unname(c(bayes_coef[2], bayes_coef[3], bayes_coef[1]))
myers_means <- sapply(x_grid, Myer_harvest, 0, myers_bayes_pars)
myers_f <- function(x,h,p) Myer_harvest(x, h, p[c("r0", "theta", "K")])
head(myers_pardist)
    deviance      K      r0 theta   stdQ
170    37.69  32.42 0.32913 2.156 0.2818
171    32.99  51.72 0.20896 2.373 0.3612
172    40.26  52.83 0.22398 2.260 0.4654
173    32.80  76.54 0.12968 2.749 0.3711
174    32.81 216.04 0.04136 3.512 0.3907
175    32.55 330.07 0.02654 3.813 0.3118
myers_bayes_pars
[1] 31.8192  0.1065 34.4760

### Phase-space diagram of the expected dynamics

models <- data.frame(x=x_grid,
GP=tgp_dat$y, True=true_means, MLE=est_means, Ricker=ricker_means, Allen = allen_means, Myers = myers_means) models <- melt(models, id="x") # some labels names(models) <- c("x", "method", "value") # labels for the colorkey too model_names = c("GP", "True", "MLE", "Ricker", "Allen", "Myers") colorkey=cbPalette names(colorkey) = model_names  plot_gp <- ggplot(tgp_dat) + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") + geom_line(data=models, aes(x, value, col=method), lwd=1, alpha=0.8) + geom_point(data=obs, aes(x,y), alpha=0.8) + xlab(expression(X[t])) + ylab(expression(X[t+1])) + scale_colour_manual(values=cbPalette) print(plot_gp) ## Goodness of fit This shows only the mean predictions. For the Bayesian cases, we can instead loop over the posteriors of the parameters (or samples from the GP posterior) to get the distribution of such curves in each case. require(MASS) step_ahead <- function(x, f, p){ h = 0 x_predict <- sapply(x, f, h, p) n <- length(x_predict) - 1 y <- c(x[1], x_predict[1:n]) y } step_ahead_posteriors <- function(x){ gp_f_at_obs <- gp_predict(gp, x, burnin=1e4, thin=300) df_post <- melt(lapply(sample(100), function(i){ data.frame(time = 1:length(x), stock = x, GP = mvrnorm(1, gp_f_at_obs$Ef_posterior[,i], gp_f_at_obs$Cf_posterior[[i]]), True = step_ahead(x,f,p), MLE = step_ahead(x,f,est$p),
}), id=c("time", "stock"))
}

ggplot(df_post) + geom_point(aes(time, stock)) +
geom_line(aes(time, value, col=variable, group=interaction(L1,variable)), alpha=.1) +
scale_colour_manual(values=colorkey, guide = guide_legend(override.aes = list(alpha = 1))) 

## Optimal policies by value iteration

Compute the optimal policy under each model using stochastic dynamic programming. We begin with the policy based on the GP model,

MaxT = 1000
# uses expected values from GP, instead of integrating over posterior
#matrices_gp <- gp_transition_matrix(gp_dat$E_Ef, gp_dat$E_Vf, x_grid, h_grid)

# Integrate over posteriors
matrices_gp <- gp_transition_matrix(gp_dat$Ef_posterior, gp_dat$Vf_posterior, x_grid, h_grid)

# Solve the SDP using the GP-derived transition matrix
opt_gp <- value_iteration(matrices_gp, x_grid, h_grid, MaxT, xT, profit, delta, reward)

Determine the optimal policy based on the allen and MLE models

matrices_true <- f_transition_matrix(f, p, x_grid, h_grid, sigma_g)
opt_true <- value_iteration(matrices_true, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)

matrices_estimated <- f_transition_matrix(est$f, est$p, x_grid, h_grid, est$sigma_g) opt_estimated <- value_iteration(matrices_estimated, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta) Determine the optimal policy based on Bayesian Allen model matrices_allen <- parameter_uncertainty_SDP(allen_f, x_grid, h_grid, pardist, 4) opt_allen <- value_iteration(matrices_allen, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta) Bayesian Ricker matrices_ricker <- parameter_uncertainty_SDP(ricker_f, x_grid, h_grid, as.matrix(ricker_pardist), 3) opt_ricker <- value_iteration(matrices_ricker, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta) Bayesian Myers model matrices_myers <- parameter_uncertainty_SDP(myers_f, x_grid, h_grid, as.matrix(myers_pardist), 4) myers_alt <- value_iteration(matrices_myers, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta) Assemble the data OPT = data.frame(GP = opt_gp$D, True = opt_true$D, MLE = opt_estimated$D, Ricker = opt_ricker$D, Allen = opt_allen$D, Myers = myers_alt$D) colorkey=cbPalette names(colorkey) = names(OPT)  ## Graph of the optimal policies policies <- melt(data.frame(stock=x_grid, sapply(OPT, function(x) x_grid[x])), id="stock") names(policies) <- c("stock", "method", "value") ggplot(policies, aes(stock, stock - value, color=method)) + geom_line(lwd=1.2, alpha=0.8) + xlab("stock size") + ylab("escapement") + scale_colour_manual(values=colorkey) ## Simulate 100 realizations managed under each of the policies sims <- lapply(OPT, function(D){ set.seed(1) lapply(1:100, function(i) ForwardSimulate(f, p, x_grid, h_grid, x0, D, z_g, profit=profit, OptTime=OptTime) ) }) dat <- melt(sims, id=names(sims[[1]][[1]])) dt <- data.table(dat) setnames(dt, c("L1", "L2"), c("method", "reps")) # Legend in original ordering please, not alphabetical: dt$method = factor(dt$method, ordered=TRUE, levels=names(OPT)) ggplot(dt) + geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=.1) + scale_colour_manual(values=colorkey, guide = guide_legend(override.aes = list(alpha = 1))) Profit <- dt[, sum(profit), by=c("reps", "method")] Profit[, mean(V1), by="method"]  method V1 1: GP 24.908 2: True 26.532 3: MLE 4.420 4: Ricker 8.363 5: Allen 7.041 6: Myers 7.347 ggplot(Profit, aes(V1)) + geom_histogram() + facet_wrap(~method, scales = "free_y") + guides(legend.position = "none") + xlab("Total profit by replicate") 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))

c(allen = allen_deviance, ricker=ricker_deviance, myers=myers_deviance, true=true_deviance, mle=mle_deviance)
  allen  ricker   myers    true     mle
45.26   44.78   34.48  -61.08 -287.60