text:
## ----plot-options, message=FALSE, warning=FALSE, include=FALSE, echo=FALSE----
cbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
               "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

## ----libraries, cache=FALSE, message=FALSE, warning=FALSE, include=FALSE, echo=FALSE----
library("methods")
library("knitr")

library("nonparametricbayes")
library("pdgControl")

library("ggplot2")
library("reshape2")
library("plyr")
library("data.table")

library("R2jags")
library("emdbook") # for as.mcmc.bugs (?)
library("coda")  # for as.mcmc

library("modeest")
library("MASS")
library("pander")

opts_chunk$set(fig.width=4, fig.height=3, echo=FALSE)
theme_set(theme_bw(base_size=12))
toggle = "hide" # results argument

## ----posterior-mode------------------------------------------------------
library(modeest)
posterior.mode <- function(x) {
  mlv(x, method="shorth")$M
}

## ----stateeq-------------------------------------------------------------
f <- RickerAllee
p <- c(2, 8, 2)
K <- 8
allee <- 2

## ----sdp-pars, dependson="stateeq"---------------------------------------
sigma_g <- 0.05
sigma_m <- 0.0
z_g <- function() rlnorm(1, 0, sigma_g)
z_m <- function() 1
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 <- 10 # simulation under policy starts from
Tobs <- 40
MaxT = 1000 # timeout for value iteration convergence

## ----obs, dependson="sdp-pars", fig.keep='none'--------------------------
  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()

## ----mle, dependson="obs"------------------------------------------------
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,.1),
         p[2]*rlnorm(1,0,.1),
         p[3]*rlnorm(1,0, .1),
         sigma_g * rlnorm(1,0,.1))
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)

## ----mle-output, dependson="mle", results=toggle-------------------------
true_means <- sapply(x_grid, f, 0, p)
est_means <- sapply(x_grid, est$f, 0, est$p)

## ----gp-priors-----------------------------------------------------------
s2.p <- c(5,5)
d.p = c(10, 1/0.1)

## ----gp, dependson=c("gp-priors", "obs")---------------------------------
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)

## ----gp_traces_densities, dependson="gp"---------------------------------
gp_assessment_plots <- summary_gp_mcmc(gp, burnin=1e4, thin=300)

## ----gp-output, dependson="gp", results=toggle---------------------------
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) )

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

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

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

## ----allen-priors, dependson="common-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)

## ----allen_priors_xtable, include=TRUE, echo=FALSE-----------------------
allen_priors_xtable <- data.frame(parameter = c("$r$", "$K$", "$X_C$", "$\\sigma$"),
           "lower bound" = c(r0_prior_p[1], K_prior_p[1], theta_prior_p[1], stdQ_prior_p[1]),
           "upper bound" = c(r0_prior_p[2], K_prior_p[2], theta_prior_p[2], stdQ_prior_p[2]))

## ----allen-mcmc, dependson=c("allen-model", "allen-pars", "jags-setup"), results=toggle----
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))

## ----allen-traces, dependson="allen-mcmc"--------------------------------
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")
plot_allen_traces <- ggplot(allen_posteriors) + geom_line(aes(index, value)) +
  facet_wrap(~ variable, scale="free", ncol=1)

## ----allen-posteriors, dependson=c("allen-traces", "allen-priors")-------
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))
})
plot_allen_posteriors <- 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)

## ----allen-output, dependson=c("posterior-mode", "allen-traces"), results=toggle----
A <- allen_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
pardist <- acast(A, index ~ variable)
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
head(pardist)

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

## ----ricker-priors, dependson="common-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)

## ----ricker_priors_xtable, include=TRUE, echo=FALSE, cache=FALSE---------
ricker_priors_xtable <- data.frame(
  parameter = c("$r$", "$K$", "$\\sigma$"),
  "lower bound" = c(r0_prior_p[1], K_prior_p[1], stdQ_prior_p[1]),
  "upper bound" = c(r0_prior_p[2], K_prior_p[2], stdQ_prior_p[2]))

## ----ricker-mcmc, dependson="ricker-model", results=toggle---------------
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"))

## ----ricker-traces, dependson="ricker-mcmc"------------------------------
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")
plot_ricker_traces <- ggplot(ricker_posteriors) + geom_line(aes(index, value)) +
  facet_wrap(~ variable, scale="free", ncol=1)

## ----ricker-posteriors, dependson=c("ricker-traces", "ricker-priors")----
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
plot_ricker_posteriors <- 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") +  # don't plot priors
  facet_wrap(~ variable, scale="free", ncol=2)

## ----ricker-output, dependson=c("posterior-mode", "ricker-traces"), results=toggle----
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)
bayes_coef <- apply(ricker_pardist,2, posterior.mode)
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)
ricker_bayes_pars

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

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

## ----myers_prior_xtable, include=TRUE, echo=FALSE------------------------
myers_priors_xtable <- data.frame(parameter = c("$r$", "$K$", "$\\theta$", "$\\sigma$"),
           "lower bound" = c(r0_prior_p[1], K_prior_p[1], theta_prior_p[1], stdQ_prior_p[1]),
           "upper bound" = c(r0_prior_p[2], K_prior_p[2], theta_prior_p[2], stdQ_prior_p[2]))

## ----myers-mcmc, dependson="myers-model", results=toggle-----------------
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"))

## ----myers-traces, dependson="myers-mcmc"--------------------------------
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")
plot_myers_traces <- ggplot(myers_posteriors) +
  geom_line(aes(index, value)) + # priors, need to fix order though
  facet_wrap(~ variable, scale="free", ncol=1)

## ----myers-posteriors, dependson="myers-traces"--------------------------
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))
})
plot_myers_posteriors <- 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") +  # Whoops, these are misaligned. see table instead
  facet_wrap(~ variable, scale="free", ncol=3)

## ----myers-output, dependson=c("posterior-mode", "myers-traces"), results=toggle----
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)
bayes_coef <- apply(myers_pardist,2, posterior.mode) # much better estimates
myers_bayes_pars <- unname(c(bayes_coef["r0"], bayes_coef["theta"], bayes_coef["K"]))
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)
myers_bayes_pars

## ----assemble-models, dependson=c("myers-output", "ricker-output", "allen-output", "gp-output", "mle-output")----
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
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
}

## ----gp-opt, dependson="gp-output"---------------------------------------
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)

## ----mle-opt, dependson="mle-output"-------------------------------------
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)

## ----allen-opt, dependson="allen-output"---------------------------------
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)

## ----ricker-opt, dependson="ricker-output"-------------------------------
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)

## ----myers-opt, dependson="myers-output"---------------------------------
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-opt, dependson=c("gp-opt", "mle-opt", "allen-opt", "ricker-opt", "myers-opt")----
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)

## ----sims, dependson="assemble-opt"--------------------------------------
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)
  )
})
# turn the list into a data.frame
keep.as.columns <- names(sims[[1]][[1]])
dat <- melt(sims, id=keep.as.columns)
sims_data <- data.table(dat)
setnames(sims_data, c("L1", "L2"), c("method", "reps"))
# Legend in original ordering please, not alphabetical:
sims_data$method = factor(sims_data$method, ordered=TRUE, levels=names(OPT))

## ----profits, dependson="sims", results=toggle---------------------------
Profit <- sims_data[, sum(profit), by=c("reps", "method")]
tmp <- dcast(Profit, reps ~ method)
#tmp$Allen <- tmp[,"Allen"] + rnorm(dim(tmp)[1], 0, 1) # jitter for plotting
tmp <- tmp / tmp[,"True"]
tmp <- melt(tmp[2:dim(tmp)[2]])
actual_over_optimal <-subset(tmp, variable != "True")

## ----dic_calc, dependson=c("posterior-mode", "myers-output", "ricker-output", "allen-output", "gp-output", "mle-output"), include=FALSE, echo=FALSE----
dic.dic <- function (x) sum(x$deviance) +  sum(x[[2]])
recompile(allen_jags)
allen_dic <- dic.dic(dic.samples(allen_jags$model, n.iter=1000, type="popt"))
recompile(ricker_jags)
ricker_dic <- dic.dic(dic.samples(ricker_jags$model, n.iter=1000, type="popt"))
recompile(myers_jags)
myers_dic <- dic.dic(dic.samples(myers_jags$model, n.iter=1000, type="popt"))
dictable <- data.frame(Allen = allen_dic + 2*length(bayes_pars),
                       Ricker = ricker_dic + 2*length(ricker_bayes_pars),
                       Myers = myers_dic + 2*length(myers_bayes_pars),
                       row.names = c("DIC"))

## ----figure_1, dependson=c("assemble-models", "par-fns", "plot-options"), fig.cap="Points show the training data of stock-size over time.  Curves show the posterior step-ahead predictions based on each of the estimated models. Observe that all models are fitting the data reasonably well.", fig.width=8, fig.height=6----
step_ahead_posteriors <- function(x){
  gp_f_at_obs <- gp_predict(gp, x, burnin=1e4, thin=300)
  df_post <- melt(lapply(sample(100, 30),
  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),
                Allen = step_ahead(x, allen_f, pardist[i,]),
                Ricker = step_ahead(x, ricker_f, ricker_pardist[i,]),
                Myers = step_ahead(x, myers_f, myers_pardist[i,]))
  }), id=c("time", "stock"))
}
df_post <- step_ahead_posteriors(x)
figure1b_posteriors <- ggplot(df_post) + geom_point(aes(time, stock)) +
  geom_line(aes(time, value, col=variable, group=interaction(L1,variable)), alpha=.1) +
  facet_wrap(~variable) +
  scale_colour_manual(values=colorkey, guide = guide_legend(override.aes = list(alpha = 1))) +
  theme(legend.position="none")
figure1b_posteriors

## ----deviances, dependson=c("dic_calc"), include=FALSE, echo=FALSE-------
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"))
xtable::xtable(rbind(dictable, aictable, bictable))

## ----Table1, dependson="deviances", include=TRUE, results="asis", echo=FALSE, cache=FALSE----
xtable::xtable(rbind(dictable, aictable, bictable), caption="Model selection scores for several common criteria all (wrongly) select the simplest model. As the true (Allen) model is not distinguishable from the simpler (Ricker) model in the region of the observed data, this error cannot be avoided regardless of the model choice criterion. This highlights the danger of model choice when the selected model will be used outside of the observed range of the data.")

## ----figure_2, dependson=c("assemble-models", "plot-options"), fig.cap="Graph of the inferred Gaussian process compared to the true process and maximum-likelihood estimated process.  Graph shows the expected value for the function $f$ under each model.  Two standard deviations from the estimated Gaussian process covariance with (light grey) and without (darker grey) measurement error are also shown.  The training data is also shown as black points.  The GP is conditioned on (0,0), shown as a pseudo-data point.", fig.width=8, fig.height=6----
x_grid_short <- x_grid[1:40]
gp_short <- gp_predict(gp, x_grid_short, burnin=1e4, thin=300)
models_posteriors <-
  melt(lapply(sample(100, 50),
              function(i){
    sample_gp <- mvrnorm(1,
                            gp_short$Ef_posterior[,i],
                            gp_short$Cf_posterior[[i]])
    data.frame(stock = x_grid_short,
               GP = sample_gp,
               y = sample_gp,
               ymin = sample_gp - 2 * sqrt(gp_short$E_Vf),
               ymax = sample_gp + 2 * sqrt(gp_short$E_Vf),
               True = sapply(x_grid_short,f,0, p),
               MLE = sapply(x_grid_short,f,0, est$p),
               Allen = sapply(x_grid_short, allen_f, 0, pardist[i,]),
               Ricker = sapply(x_grid_short, ricker_f, 0, ricker_pardist[i,]),
               Myers = sapply(x_grid_short, myers_f, 0, myers_pardist[i,]))
             }),
       id=c("stock", "y", "ymin", "ymax"))
ggplot(models_posteriors) +
    geom_ribbon(aes(x=stock, y=y, ymin=ymin, ymax=ymax, group=L1),
                  fill = "gray80",
                  data=subset(models_posteriors, variable == "GP")) +
    geom_line(aes(stock, value, col = variable,
                  group=interaction(L1,variable)),
              alpha=.2) +
    geom_point(data = obs, aes(x,y), alpha = 0.8) +
    xlab(expression(X[t])) + ylab(expression(X[t+1])) +
    facet_wrap(~variable) +
    scale_colour_manual(values=colorkey) +
    theme(legend.position="none")

## ----figure_3, fig.width=8, fig.height=6, fig.cap="Beyond the tipping point (outside the range of the training data), the true dynamics fall outside the uncertainty bands of the wrong parametric models, but inside the uncertainty predicted by the GP. Points show the stock size simulated by the true model.  Overlay shows the range of states predicted by each model, based on the state observed in the previous time step. The Ricker model always (wrongly) predicts positive population growth, while the actual population shrinks in each step as the initial condition falls below the Allee threshold of the underlying model (Allen).  Note that the GP is both more pessimistic and more uncertain about the future state than the parametric models.", dependson="plot-options"----
y <- numeric(8)
y[1] <- 4.5
for(t in 1:(length(y)-1))
      y[t+1] = z_g() * f(y[t], h=0, p=p)
# predicts means, does not reflect uncertainty estimate!
crash_data <- step_ahead_posteriors(y)
crash_data <- subset(crash_data, variable %in% c("GP", "Allen", "Ricker", "Myers"))
ggplot(crash_data) +
  geom_boxplot(aes(as.factor(as.integer(time)), value,
                   fill = variable, col=variable),
               alpha=.7, outlier.size=1, position="identity") +
#  geom_line(aes(time, value, col = variable,
#            group=interaction(L1,variable)), alpha=.1) +
  geom_point(aes(time, stock), size = 3) +
  scale_fill_manual(values=colorkey[c("GP", "Allen", "Ricker", "Myers")],
                      guide = guide_legend(override.aes = list(alpha = 1))) +
  scale_colour_manual(values=colorkey[c("GP", "Allen", "Ricker", "Myers")],
                      guide = guide_legend(override.aes = list(alpha = 1))) +
  facet_wrap(~variable) +
  theme(legend.position="none") + xlab("time") + ylab("stock size")

## ----figure_4, dependson = c("assemble-opt"), fig.cap="The steady-state optimal policy (infinite boundary) calculated under each model.  Policies are shown in terms of target escapement, $S_t$, as under models such as this a constant escapement policy is expected to be optimal [@Reed1979]."----
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)

## ----figure_5, dependson=c("sim"), fig.cap="Gaussian process inference outperforms parametric estimates. Shown are 100 replicate simulations of the stock dynamics (eq 1) under the policies derived from each of the estimated models, as well as the policy based on the exact underlying model.", fig.width=6, fig.height=6----
ggplot(sims_data) +
  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))) +
  facet_wrap(~method) + guides(legend.position="none")

## ----figure_6, dependson=c("profits"), fig.cap="Histograms of the realized net present value of the fishery over a range of simulated data and resulting parameter estimates. For each data set, the three models are estimated as described above. Values plotted are the averages of a given policy over 100 replicate simulations. Details and code provided in the supplement.", fig.width=6, fig.height=6----

ggplot(actual_over_optimal, aes(value)) + geom_histogram(aes(fill=variable)) +
  facet_wrap(~variable, ncol=2)  +
  guides(legend.position = "none") +
  xlab("Total profit by replicate") +
  scale_fill_manual(values=colorkey) # density plots fail when delta fn

# ggplot(actual_over_optimal, aes(value)) + geom_histogram(aes(fill=variable), binwidth=0.1) +
#  xlab("Total profit by replicate")+ scale_fill_manual(values=colorkey)
# ggplot(actual_over_optimal, aes(value, fill=variable, color=variable)) + # density plots fail when delta fn
#  stat_density(aes(y=..density..), position="stack", adjust=3, alpha=.9) +
#  xlab("Total profit by replicate")+ scale_fill_manual(values=colorkey)+ scale_color_manual(values=colorkey)

## ----figure_S1, fig.cap = "Traces from the MCMC estimates of the GP model show reasonable mixing (no trend) and sampling rejection rate (no piecewise jumps)", dependson="plot-options"----
gp_assessment_plots$traces_plot

## ----figure_S2, fig.cap="Posterior distributions from the MCMC estimate of the GP model. Prior curves shown in red."----
gp_assessment_plots$posteriors_plot

## ----figure_S3, fig.height=6, fig.cap="Traces from the MCMC estimates of the Ricker model show reasonable mixing (no trend) and sampling rejection rate (no piecewise jumps)", fig.width=6----
plot_ricker_traces

## ----figure_S4, fig.cap="Posteriors from the MCMC estimate of the Ricker model", fig.width=6, fig.height=4----
ggplot(ricker_posteriors, aes(value)) +
  stat_density(geom="path", position="identity", alpha=0.7) +
  facet_wrap(~ variable, scale="free", ncol=2)

## ----Table S1, results = "asis"------------------------------------------
pander::pandoc.table(ricker_priors_xtable,
  caption = "Parameterization range for the uniform priors in the Ricker model")

## ----figure_S5, fig.height=6, fig.cap="Traces from the MCMC estimates of the Myers model show reasonable mixing (no trend) and sampling rejection rate (no piecewise jumps)", fig.width=6----
plot_myers_traces

## ----figure_S6, fig.cap="Posterior distributions from the MCMC estimates of the Myers model", fig.width=6, fig.height=6----
ggplot(myers_posteriors, aes(value)) +
  stat_density(geom="path", position="identity", alpha=0.7) +
  facet_wrap(~ variable, scale="free", ncol=2)

## ----TableS2, results="asis"---------------------------------------------
pander::pandoc.table(myers_priors_xtable,
           caption = "Parameterization range for the uniform priors in the Myers model")

## ----figure_S7, fig.height=6, fig.cap="Traces from the MCMC estimates of the Allen model show reasonable mixing (no trend) and sampling rejection rate (no piecewise jumps)", fig.width=6----
plot_allen_traces

## ----figure_S8, fig.cap="Posteriors from the MCMC estimate of the Allen model", fig.width=6, fig.height=6----
ggplot(allen_posteriors, aes(value)) +
  stat_density(geom="path", position="identity", alpha=0.7) +
  facet_wrap(~ variable, scale="free", ncol=2)

## ----TableS3, results = "asis"-------------------------------------------
pander::pandoc.table(allen_priors_xtable,
  caption = "Parameterization range for the uniform priors in the Allen model")

## ----include=FALSE-------------------------------------------------------
unlink("ricker_process.bugs")
unlink("allen_process.bugs")
unlink("myers_process.bugs")

save(list=ls(), file="small-allee.rda")


text:
```{r plotting-options, message=FALSE, warning=FALSE, include=FALSE}
opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE, cache=FALSE, comment=NA)
library(ggplot2) # plotting
opts_knit$set(upload.fun =imgur_upload)
theme_set(theme_bw(base_size=10))
theme_update(panel.background = element_rect(fill = "transparent", colour = NA),
             plot.background = element_rect(fill = "transparent", colour = NA))

```

Basic regression in Gaussian processes
---------------------------------------

Working out coded examples for basic Gaussian process regression using R.  I've just read through the first few chapters of [Rasmussen & Williams (2006)](http://www.GaussianProcess.org/gpml), this implements the examples discussed in Chapter 2.1-2.5.

<!-- Also motivated by [James Keirstead](http://www.jameskeirstead.ca/r/gaussian-process-regression-with-r/), with cleaner code.  -->

## Function definitions

```{r}
library(mcmc)
#' Basic regression in Gaussian processes
#'
#'
#' @param x Observed x values, (vector or matrix with columns for each dimension of data)
#' @param y Vector of observed y values in the training data
#' @param init_pars the initial guesses for lengthscale l and process noise sigma_n
#' @param n iterations of the metropolis algorithm
#' @details Currently assumes the covariance function.  By default we will use
#' the squared exponential (also called radial basis or Gaussian,
#' though it is not this that gives Gaussian process it's name;
#' any covariance function would do) as the covariance function,
#' $\operatorname{cov}(X_i, X_j) = \exp\left(-\frac{(X_i - X_j)^2}{2 \ell ^ 2}\right)$
#' @return the MCMC output, as constructed by \link{metrop} from the mcmc package.
#' @export
#' @import mcmc
gp_mcmc <- function(x, y, init_pars = c(l=1, sigma.n=1), n = 1e4, d.p = c(5,5), s2.p = c(5,5)){

  lpriors <- function(pars){
    prior <- unname(
      dgamma(exp(pars[1]), d.p[1], scale = d.p[2]) *
      dgamma(exp(pars[2]), s2.p[1], s2.p[2])
    )
    log(prior)
  }
  posterior <- function(pars, x, y){
    l <- exp(pars[1])
    sigma.n <- exp(pars[2])

    cov <- function(X, Y) outer(X, Y, SE, l)
    I <- diag(1, length(x))
    K <- cov(x, x)

    loglik <- - 0.5 * t(y) %*% solve(K + sigma.n^2 * I) %*% y -
      log(det(K + sigma.n^2*I)) -
      length(y) * log(2 * pi) / 2

    loglik + lpriors(pars)
  }

  out <- metrop(posterior, log(init_pars), n, x = obs$x, y = obs$y)
  out$d.p <- d.p
  out$s2.p <- s2.p

  out
}
```

```{r}
#' predict the expected values and posterior distributions of the Gaussian Process
#'
#' @param x_predict the values at which we desire predictions
#' @param covs if TRUE, return covariances instead of variances (set high thinning as this is memory intensive)
#' @import reshape2
#' @export
gp_predict <- function(gp, x_predict, burnin=0, thin=1, covs=FALSE){

  postdist <- cbind(index=1:gp$nbatch, as.data.frame(exp(gp$batch)))
  s <- seq(burnin+1, gp$nbatch, by=thin)
  postdist <- postdist[s,]
  names(postdist) <- c("index", names(gp$initial))

  indices <- 1:dim(postdist)[1]

  out <- lapply(indices, function(i){
    l <- postdist[i, "l"]
    sigma.n <- postdist[i, "sigma.n"]

    cov <- function(X, Y) outer(X, Y, SE, l)
    cov_xx_inv <- solve(cov(obs$x, obs$x) + sigma.n^2 * diag(1, length(obs$x)))
    Ef <- cov(x_predict, obs$x) %*% cov_xx_inv %*% obs$y
    Cf <- cov(x_predict, x_predict) - cov(x_predict, obs$x)  %*% cov_xx_inv %*% cov(obs$x, x_predict)
    list(Ef = Ef, Cf = Cf, Vf = diag(Cf))
  })

  Ef_posterior <- sapply(out, `[[`, "Ef")
  Cf_posterior <- sapply(out, `[[`, "Cf")
  Vf_posterior <- sapply(out, `[[`, "Vf")

  E_Ef <- rowMeans(Ef_posterior)
  E_Cf <- matrix( apply(Cf_posterior, 1, sum) / dim(Cf_posterior)[2], ncol = sqrt(dim(Cf_posterior)[1]) )
  E_Vf <- diag(E_Cf) # same as rowMeans(Vf_posterior)

  # list format better for return
  Cf_posterior <- lapply(out, `[[`, "Cf")

  list(Ef_posterior = Ef_posterior,
       Vf_posterior = Vf_posterior,
       Cf_posterior = Cf_posterior,
       E_Ef = E_Ef, E_Cf = E_Cf, E_Vf = E_Vf)

}
```

```{r}
library(reshape2)
library(ggplot2)
library(plyr)
#' Summary plots showing the trace and posteriors for the gp_mcmc estimates
#'
#' @param gp a fit of the gaussian process from gp_mcmc
#' @param burnin length of sequence to discard as a transient
#' @param thin frequency of sub-sampling (make posterior distribution smaller if necessary)
#' @return two ggplot2 objects, one plotting the trace and one
#' plotting the posteriors in black with priors overlaid in red.
#' @import reshape2
#' @import ggplot2
#' @import plyr
#' @export
summary_gp_mcmc <- function(gp, burnin=0, thin=1){

  postdist <- cbind(index=1:gp$nbatch, as.data.frame(exp(gp$batch)))
  s <- seq(burnin+1, gp$nbatch, by=thin)
  postdist <- postdist[s,]
  names(postdist) <- c("index", names(gp$initial))

  # TRACES
  df <- melt(postdist, id="index")
  traces_plot <-
    ggplot(df) + geom_line(aes(index, value)) +
    facet_wrap(~ variable, scale="free", ncol=1)

  s2.p <- gp$s2.p
  d.p <- gp$d.p
  s2_prior <- function(x) dgamma(x, s2.p[1], s2.p[2])
  d_prior <- function(x) dgamma(x, d.p[1], scale = d.p[2])
  prior_fns <- list(l = d_prior, sigma.n = s2_prior)

  prior_curves <- ddply(df, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = prior_fns[[dd$variable[1]]](grid))
  })

  # Posteriors (easier to read than histograms)
  posteriors_plot <- ggplot(df, aes(value)) +
    stat_density(geom="path", position="identity", alpha=0.7) +
    geom_line(data=prior_curves, aes(x=value, y=density), col="red") +
    facet_wrap(~ variable, scale="free", ncol=2)

  print(traces_plot)
  print(posteriors_plot)

  out <- list(traces_plot = traces_plot, posteriors_plot = posteriors_plot)
  invisible(out)
}

# Helper function: The default covariance function
SE <- function(Xi,Xj, l) exp(-0.5 * (Xi - Xj) ^ 2 / l ^ 2)
```

## Example

```{r}
x_predict <- seq(-5,5,len=50)
obs <- data.frame(x = c(-4, -3, -1,  0,  2),
                  y = c(-2,  0,  1,  2, -1))
```

```{r}
gp <- gp_mcmc(obs$x, obs$y, d.p = c(2,.5), s2.p = c(2,.5))
```

```{r}
summary_gp_mcmc(gp)
```

```{r}
gp_dat <- gp_predict(gp, x_predict)
```

```{r}
tgp_dat <-
    data.frame(  x = x_predict,
                 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) )

ggplot(tgp_dat) + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") +
    geom_line(aes(x, y), lwd=2, alpha=0.8) +
    geom_point(data=obs, aes(x,y), alpha=0.8) +
    xlab(expression(X[t])) + ylab(expression(X[t+1]))
```

repo_name: cboettig/nonparametric-bayes

repo_language: TeX

mime_type: text/plain


text:
Generate data under Beverton Holt model
==============================================

```{r libraries}
require(pdgControl)
require(nonparametricbayes)
require(ggplot2)
```

### Beverton-Holt function

Simulate some training data under a stochastic growth function with standard parameterization,

``` {r stateeqn}
f <- BevHolt
p <- c(1.5,.05)
K <- (p[1]-1)/p[2]
```

Parameter definitions

```{r}
x_grid = seq(0, 1.5 * K, length=101)
T <- 40
sigma_g <- 0.1
x <- numeric(T)
x[1] <- 1
```

Noise function, profit function
```{r}
z_g <- function() rlnorm(1, 0, sigma_g) #1+(2*runif(1, 0,  1)-1)*sigma_g #
profit <- profit_harvest(1,0,0)
```

Simulation

```{r}
set.seed(111)
for(t in 1:(T-1))
  x[t+1] = z_g() * f(x[t], h=0, p=p)
```

Predict the function over the target grid (lag-1)

```{r}
obs <- data.frame(x=x[1:(T-1)],y=x[2:T])
X <- x_grid
```


text:
```{r plotting-options}
opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE, cache=FALSE, comment=NA, dev="CairoPDF", fig.width=6, fig.height=4)
library(ggplot2) # plotting
opts_knit$set(upload.fun = socialR::flickr.url)
theme_set(theme_bw(base_size=10))
theme_update(panel.background = element_rect(fill = "transparent", colour = NA),
             plot.background = element_rect(fill = "transparent", colour = NA))
cbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
               "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
```

```{r}
set.seed(12345)
```

```{r}

#' @import MASS
#' @import reshape2
samples <- function(obs, l=1, sigma.n=0.3, x_predict = seq(-5,5,len=50), n=3){

  require(MASS); require(reshape2)
  SE <- function(Xi,Xj, l) exp(-0.5 * (Xi - Xj) ^ 2 / l ^ 2)
  cov <- function(X, Y) outer(X, Y, SE, l)
  cov_xx_inv <- solve(cov(obs$x, obs$x) + sigma.n^2 * diag(1, length(obs$x)))
  Ef <- cov(x_predict, obs$x) %*% cov_xx_inv %*% obs[[2]]
  Cf <- cov(x_predict, x_predict) -
        cov(x_predict, obs$x)  %*% cov_xx_inv %*% cov(obs$x, x_predict)

  values <- mvrnorm(n, Ef, Cf)
  dat <- data.frame(x=x_predict, t(values))
  dat <- melt(dat, id="x")
  Ey <- data.frame(x=x_predict,
                   y=(Ef),
                   ymin=(Ef-2*sqrt(diag(Cf))),
                   ymax=(Ef+2*sqrt(diag(Cf))))
  list(dat=dat, Ey=Ey, obs=obs, x_predict=x_predict, l=l, sigma.n=sigma.n, n=n)

  merge(obs, merge(dat, Ey, "x", all=TRUE), "x", all=TRUE)
}

x_predict <- seq(-5,5,len=50)
obs <- data.frame(x = c(x_predict[c(5,10,15,20,25)]),
                  observed = x_predict[c(10,  15,  20,  25, 10)])

S <- samples(obs, l=1, sigma.n=.5, x_predict, n=3)

ggplot(S) +
  geom_ribbon(aes(x,y,ymin=ymin, ymax=ymax), fill="grey80") +
  geom_point(aes(x,observed)) +
  geom_line(aes(x,value, color=variable))
```

```{r}
ggplot(S$dat,aes(x=x,y=value)) +
  geom_ribbon(data=S$Ey, aes(x=x,y=y, ymin=ymin, ymax=ymax), fill="grey80") + # Var
  geom_line(aes(color=variable)) + #REPLICATES
  geom_line(data=S$Ey,aes(x=x,y=y), size=1) + #MEAN
  geom_point(data=S$obs,aes(x=x,y=y)) +  #OBSERVED DATA
  scale_y_continuous(lim=c(-4,4), name="output, f(x)") +
  xlab("input, x")
```

Note that unlike the previous case, the posterior no longer collapses completely around the neighborhood of the test points.

We can also compute the likelihood (and marginal likelihood over the prior) of the data directly from the inferred multivariate normal distribution, which can allow us to tune the hyperparameters such as the characteristic length scale $\ell$ and the observation noise $\sigma_n$.  The most obvious approach would be to do so by maximum likelihood, giving point estimates of the hyper-parameters, though presumably we could be Bayesian about these as well.

<!-- Continuing... 2012-10-22 -->

The likelihood is given by
$$\log(p(y | X)) = -\tfrac{1}{2} \mathbf{y}^T (K + \sigma_n^2 \mathbf{I})^{-1} y  - \tfrac{1}{2} \log\left| K + \sigma_n^2 \mathbf{I} \right| - \tfrac{n}{2}\log 2 \pi$$

Which is:

```{r}

#' observation pairs, obs$x, obs$x,
#' SE covariance between two observations
minusloglik <- function(pars){

  l <- pars[1]
  sigma.n <- pars[2]

  cov <- function(X, Y) outer(X, Y, SE, l)
  I <- diag(1, length(obs$x))
  K <- cov(obs$x, obs$x)

  0.5 * t(obs$y) %*% solve(K + sigma.n^2 * I) %*% obs$y +
    log(det(K + sigma.n^2*I)) +
    length(obs$y) * log(2 * pi) / 2
  }
```

```{r}
pars <- c(l=1, sigma.n=1)
o <- optim(pars, minusloglik)
o$par
```

```{r}
lpriors <- function(pars){
   d.p <- c(5, 5)
  s2.p <- c(5, 5)

  prior <- unname(
    dgamma(exp(pars[1]), d.p[1], scale = d.p[2]) *
    dgamma(exp(pars[2]), s2.p[1], s2.p[2])
#    dunif(exp(pars[1]), 0, 100) *
#    dunif(exp(pars[2]), 0, 100)
  )

  log(prior)
}

posterior <- function(pars, x, y){

  l <- exp(pars[1])
  sigma.n <- exp(pars[2])

  cov <- function(X, Y) outer(X, Y, SE, l)
  I <- diag(1, length(x))
  K <- cov(x, x)

  loglik <- - 0.5 * t(y) %*% solve(K + sigma.n^2 * I) %*% y -
    log(det(K + sigma.n^2*I)) -
    length(y) * log(2 * pi) / 2

  loglik + lpriors(pars)
}

```

```{r}
posterior(log(pars), obs$x, obs$y)
posterior(log(o$par), obs$x, obs$y)
```

```{r}
require(mcmc)
n <- 1e4
out <- metrop(posterior, log(pars), n, x = obs$x, y = obs$y)
out$accept
```

```{r gp_mcmc}
postdist <- cbind(index=1:n, as.data.frame(exp(out$batch)))
names(postdist) <- c("index", names(pars))
df <- melt(postdist, id="index")
# TRACES
ggplot(df) + geom_line(aes(index, value)) + facet_wrap(~ variable, scale="free", ncol=1)

ggplot(df) + geom_line(aes(index, log(value))) + facet_wrap(~ variable, scale="free", ncol=1)
```

```{r gp_posteriors}

d.p <- c(5, 5)
s2.p <- c(5, 5)

s2_prior <- function(x) dgamma(x, s2.p[1], s2.p[2])
d_prior <- function(x) dgamma(x, d.p[1], scale = d.p[2])
prior_fns <- list(l = d_prior, sigma.n = s2_prior)

require(plyr)
prior_curves <- ddply(df, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = prior_fns[[dd$variable[1]]](grid))
})

# Posteriors (easier to read than histograms)
ggplot(df, aes(value)) +
  stat_density(geom="path", position="identity", alpha=0.7) +
  geom_line(data=prior_curves, aes(x=value, y=density), col="red") +
  facet_wrap(~ variable, scale="free", ncol=2)
```


text:
GP Example using the May (1979) bistable model
==============================================

```
## Loading required package: bibtex
```

```
## Warning: replacing previous import 'write.bib' when loading 'pkgmaker'
```

```r
f <- May
p <- c(r = .75, k = 10, a=1.3, H=1, Q = 3)
K <- 8
```

We use the model of May _et. al._ (1979).

```r
sigma_g <- 0.04
z_g <- function(sigma_g) rlnorm(1, 0, sigma_g) #1+(2*runif(1, 0,  1)-1)*sigma_g #
x_grid <- seq(0, 1.5 * K, length=101)
h_grid <- x_grid
profit = function(x,h) pmin(x, h)
delta <- 0.01
OptTime = 20
reward = profit(x_grid[length(x_grid)], x_grid[length(x_grid)]) + 1 / (1 - delta) ^ OptTime
```

With parameters `0.75, 10, 1.3, 1, 3`.

```r
xT <- x_grid[2]
x_0_observed <- x_grid[60]
```

```r
Tobs <- 100
x <- numeric(Tobs)
x[1] <- x_0_observed
for(t in 1:(Tobs-1))
  x[t+1] = z_g(sigma_g) * f(x[t], h=0, p=p)
plot(x)
```

![plot of chunk sim-obs](http://carlboettiger.info/assets/figures/2012-12-15-a3aad15b15-sim-obs.png)

We simulate data under this model, starting from a size of `7.08`.

```r
obs <- data.frame(x=c(0,x[1:(Tobs-1)]),y=c(0,x[2:Tobs]))
```

We consider the observations as ordered pairs of observations of current stock size $x_t$ and observed stock in the following year, $x_{t+1}$.  We add the pseudo-observation of $0,0$.  Alternatively we could condition strictly on solutions passing through the origin, though in practice the weaker assumption is often sufficient.

```r
estf <- function(p){
  mu <- log(obs$x) + p["r"]*(1-obs$x/p["K"])
  -sum(dlnorm(obs$y, mu, p["s"]), log=TRUE)
}
o <- optim(par = c(r=1,K=mean(x),s=1), estf, method="L", lower=c(1e-3,1e-3,1e-3))
f_alt <- Ricker
p_alt <- c(o$par['r'], o$par['K'])
sigma_g_alt <- o$par['s']
```

```r
gp <- bgp(X=obs$x, XX=x_grid, Z=obs$y, verb=0,
          meanfn="linear", bprior="b0", BTE=c(2000,6000,2), m0r1=FALSE,
          corr="exp", trace=TRUE, beta = c(0,0),
          s2.p = c(50,50), d.p = c(10, 1/0.01, 10, 1/0.01), nug.p = c(10, 1/0.01, 10, 1/0.01),
          s2.lam = "fixed", d.lam = "fixed", nug.lam = "fixed",
          tau2.lam = "fixed", tau2.p = c(50,1))
```

We fit a Gaussian process with

```r
V <- gp$ZZ.ks2
Ef = gp$ZZ.km
tgp_dat <- data.frame(x   = gp$XX[[1]],
                  y   = gp$ZZ.km,
                 ymin = gp$ZZ.km - 1.96 * sqrt(gp$ZZ.ks2),
                 ymax = gp$ZZ.km + 1.96 * sqrt(gp$ZZ.ks2))
```

```r
true <- data.frame(x=x_grid, y=sapply(x_grid,f, 0, p))
ggplot(tgp_dat)  + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") +
  geom_line(aes(x,y)) + geom_point(data=obs, aes(x,y)) +
  geom_line(data=true, aes(x,y), col='red', lty=2)
```

![plot of chunk gp-plot](http://carlboettiger.info/assets/figures/2012-12-15-a3aad15b15-gp-plot.png)

## The transition matrix of the inferred process

```r
X <- numeric(length(x_grid))
X[38] = 1
h <- 0
```

```r
F_ <- gp_F(h, Ef, V, x_grid)
xt1 <- X %*% F_
xt10 <- xt1
for(s in 1:OptTime)
  xt10 <- xt10 %*% F_
qplot(x_grid, xt10[1,]) + geom_point(aes(y=xt1[1,]), col="grey")
```

![plot of chunk gp-F-sim](http://carlboettiger.info/assets/figures/2012-12-15-a3aad15b15-gp-F-sim.png)

```r
F_true <- par_F(h, f, p, x_grid, sigma_g)
yt1 <- X %*% F_true
yt10 <- yt1
for(s in 1:OptTime)
  yt10 <- yt10 %*% F_true
qplot(x_grid, yt10[1,]) + geom_point(aes(y=yt1[1,]), col="grey")
```

![plot of chunk par-F-sim](http://carlboettiger.info/assets/figures/2012-12-15-a3aad15b15-par-F-sim.png)

```r
transition <- melt(data.frame(x = x_grid, gp = xt1[1,], parametric = yt1[1,]), id="x")
ggplot(transition) + geom_point(aes(x,value, col=variable))
```

![plot of chunk F-sim-plot](http://carlboettiger.info/assets/figures/2012-12-15-a3aad15b15-F-sim-plot.png)

```r
matrices_gp <- gp_transition_matrix(Ef, V, x_grid, h_grid)
opt_gp <- find_dp_optim(matrices_gp, x_grid, h_grid, OptTime, xT, profit, delta, reward=reward)
```

```r
matrices_true <- f_transition_matrix(f, p, x_grid, h_grid, sigma_g)
opt_true <- find_dp_optim(matrices_true, x_grid, h_grid, OptTime, xT, profit, delta=delta, reward = reward)
```

```r
matrices_estimated <- f_transition_matrix(f_alt, p_alt, x_grid, h_grid, sigma_g_alt)
opt_estimated <- find_dp_optim(matrices_estimated, x_grid, h_grid, OptTime, xT, profit, delta=delta, reward = reward)
```

```r
policies <- melt(data.frame(stock=x_grid,
                            GP = x_grid[opt_gp$D[,1]],
                            Exact = x_grid[opt_true$D[,1]],
                            Approx = x_grid[opt_estimated$D[,1]]),
                  id="stock")

policy_plot <- ggplot(policies, aes(stock, stock - value, color=variable)) +
  geom_point() + xlab("stock size") + ylab("escapement")
policy_plot
```

![plot of chunk policy_plot](http://carlboettiger.info/assets/figures/2012-12-15-a3aad15b15-policy_plot.png)

```r
z_g <- function() rlnorm(1,0, sigma_g)
```

```r
set.seed(1)
sim_gp <- ForwardSimulate(f, p, x_grid, h_grid, K, opt_gp$D, z_g, profit=profit)
set.seed(1)
sim_true <- ForwardSimulate(f, p, x_grid, h_grid, K, opt_true$D, z_g, profit=profit)
set.seed(1)
sim_est <- ForwardSimulate(f, p, x_grid, h_grid, K, opt_estimated$D, z_g, profit=profit)
```

```r
dat <- list(est = sim_est, gp = sim_gp, true = sim_true)
dat <- melt(dat, id=names(dat[[1]]))
dt <- data.table(dat)
setnames(dt, "L1", "method")
```

```r
ggplot(dt) + geom_line(aes(time,fishstock, color=method))
```

![plot of chunk sim-fish](http://carlboettiger.info/assets/figures/2012-12-15-a3aad15b15-sim-fish.png)

```r
ggplot(dt) + geom_line(aes(time,harvest, color=method))
```

![plot of chunk sim-harvest](http://carlboettiger.info/assets/figures/2012-12-15-a3aad15b15-sim-harvest.png)

```r
c( gp = sum(sim_gp$profit), true = sum(sim_true$profit), est = sum(sim_est$profit))
```

```
    gp   true    est
 8.469 19.680 10.200
```

May RM, Beddington JR, Clark CW, Holt SJ and Laws RM (1979).
“Management of Multispecies Fisheries.”
Science, 205.
ISSN 0036-8075, http://dx.doi.org/10.1126/science.205.4403.267.


text:
## Optimization over hyperparameters:

```r
opts_knit$set(upload.fun = socialR::notebook.url)
```

Static definition of gp function

```r
gp_fit <- function(obs, X, pars = c(sigma_n = 1, l = 1)) {
    sigma_n <- pars["sigma_n"]
    l <- pars["l"]

    SE <- function(Xi, Xj, l = l) exp(-0.5 * (Xi - Xj)^2/l^2)
    cov <- function(X, Y) outer(X, Y, SE, l)

    n <- length(obs$x)
    K <- cov(obs$x, obs$x)
    I <- diag(1, n)

    cov_xx_inv <- solve(K + sigma_n^2 * I)
    Ef <- cov(X, obs$x) %*% cov_xx_inv %*% obs$y
    Cf <- cov(X, X) - cov(X, obs$x) %*% cov_xx_inv %*% cov(obs$x, X)
    llik <- -0.5 * t(obs$y) %*% cov_xx_inv %*% obs$y - 0.5 * log(det(cov_xx_inv)) -
        n * log(2 * pi)/2

    out <- list(Ef = Ef, Cf = Cf, llik = llik, obs = obs, X = X, pars = pars,
        K = K, inv = cov_xx_inv)
    class(out) = "gpfit"
    out
}

plot.gpfit <- function(gp) {
    dat <- data.frame(x = gp$X, y = (gp$Ef), ymin = (gp$Ef - 2 * sqrt(diag(gp$Cf))),
        ymax = (gp$Ef + 2 * sqrt(diag(gp$Cf))))
    ggplot(dat) + geom_ribbon(aes(x = x, y = y, ymin = ymin, ymax = ymax), fill = "grey80") +
        geom_line(aes(x = x, y = y), size = 1) + geom_point(data = gp$obs, aes(x = x,
        y = y)) + labs(title = paste("llik =", prettyNum(gp$llik)))
}
```

Estimate a GP from data with given hyperparameters, then redraw data from an actual multivariate normal.  Presumably optimal hyperparameters might be near the given values then.

```r
require(MASS)
X <- seq(-5, 5, len = 20)
obs <- data.frame(x = c(-4, -3, -1, 0, 2, 3), y = c(-2, 0, 1, 2,
    -1, -1))

fit <- gp_fit(obs, X, pars = c(sigma_n = 0.3, l = 0.6))
gp_dat <- mvrnorm(1, fit$Ef, fit$Cf)
obs <- data.frame(x = X, y = gp_dat)
X <- seq(-5, 5, len = 50)
```

Example plot

```r
plot.gpfit(gp_fit(obs, X, pars = c(sigma_n = 0.3, l = 0.6)))
```

![plot of chunk unnamed-chunk-4](http://carlboettiger.info/assets/figures/2012-11-15-28b3256dfe-unnamed-chunk-4.png)

estimates non-boundary solution of `l`:

```r
minusloglik <- function(par) {
    gp <- gp_fit(obs, X, c(sigma_n = 0.3, l = par))
    -gp$llik
}
minusloglik(1)
```

```
##       [,1]
## [1,] 37.76
```

```r
o <- optimize(minusloglik, c(0, 40))
o
```

```
## $minimum
## [1] 0.3614
##
## $objective
##       [,1]
## [1,] 26.87
```

But makes `sigma_n` arbitrarily large (boundary solution).

```r
minusloglik <- function(par) {
    gp <- gp_fit(obs, X, c(sigma_n = par, l = 0.6))
    -gp$llik
}
minusloglik(1)
```

```
##       [,1]
## [1,] 16.88
```

```r
o <- optimize(minusloglik, c(0, 40))
o
```

```
## $minimum
## [1] 40
##
## $objective
##       [,1]
## [1,] -55.4
```

```r
gp_k <- gausspr(obs$x, obs$y, kernel = "rbfdot")
```

```
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
```

```r
gp_k
```

```
## Gaussian Processes object of class "gausspr"
## Problem type: regression
##
## Gaussian Radial Basis kernel function.
##  Hyperparameter : sigma =  4.79808910472974
##
## Number of training instances learned : 20
## Train error : 0.131915781
```

Convert to `l`, `sigma=1/(2l^2)`

```r
1/sqrt(2 * as.numeric(gp_k@kernelf@kpar))
```

```
## [1] 0.3228
```

Weirdly inconsistent?

```r
gp_k <- gausspr(obs$x, obs$y, kernel = "rbfdot")
```

```
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
```

```r
1/sqrt(2 * as.numeric(gp_k@kernelf@kpar))
```

```
## [1] 0.1684
```

```r
gp_k <- gausspr(obs$x, obs$y, kernel = "rbfdot")
```

```
## Using automatic sigma estimation (sigest) for RBF or laplace kernel
```

```r
1/sqrt(2 * as.numeric(gp_k@kernelf@kpar))
```

```
## [1] 0.3206
```

### a second example, using the beverton holt data:

```r
knit("beverton_holt_data.Rmd")
```

```
##
##
## processing file: beverton_holt_data.Rmd
```

![plot of chunk unnamed-chunk-10](http://carlboettiger.info/assets/figures/2012-11-15-28b3256dfe-unnamed-chunk-10.png)

```
## output file:
## /home/cboettig/Documents/code/nonparametric-bayes/inst/examples/beverton_holt_data.md
```

```r
minusloglik <- function(par) {
    if (any(par < 0))
        Inf else {
        gp <- gp_fit(obs, X, par)
        -gp$llik
    }
}
par <- c(sigma_n = 1, l = 1)
minusloglik(par)
```

```
##       [,1]
## [1,] 184.9
```

```r
# o <- optim(par, minusloglik)
o <- optim(par, minusloglik, method = "L", lower = c(0, 0), upper = c(20,
    20))
o
```

```
## $par
## sigma_n       l
##      20      20
##
## $value
## [1] -77.64
##
## $counts
## function gradient
##        2        2
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"
```

```r
hyperpars <- o$par
```

**Yikes.**

let's try 1-D optimization:

```r
minusloglik <- function(par) {
    gp <- gp_fit(obs, X, c(sigma_n = par, l = 5.7))
    -gp$llik
}
minusloglik(1)
```

```
##      [,1]
## [1,]  105
```

```r
o <- optimize(minusloglik, c(0, 50))
o
```

```
## $minimum
## [1] 50
##
## $objective
##        [,1]
## [1,] -116.2
```

```r
hyperpars <- c(sigma_n = o$minimum, l = 1)
```

**Yikes, that shouldn't happen either**

Plot the optimal hyperparameter solution:

```r
gp <- gp_fit(obs, X, hyperpars)
df <- data.frame(x = X, y = gp$Ef, ymin = (gp$Ef - 2 * sqrt(abs(diag(gp$Cf)))),
    ymax = (gp$Ef + 2 * sqrt(abs(diag(gp$Cf)))))
true <- data.frame(x = X, y = sapply(X, f, 0, p))
ggplot(df) + geom_ribbon(aes(x, y, ymin = ymin, ymax = ymax), fill = "gray80") +
    geom_line(aes(x, y)) + geom_point(data = obs, aes(x, y)) + geom_line(data = true,
    aes(x, y), col = "red", lty = 2)
```

![plot of chunk unnamed-chunk-13](http://carlboettiger.info/assets/figures/2012-11-15-28b3256dfe-unnamed-chunk-13.png)


text:
```{r plotting-options, cache=FALSE, message=FALSE, warning=FALSE, include=FALSE}
rm(list=ls())
source("~/.knitr_defaults.R")
opts_knit$set(upload.fun = socialR::flickr.url)
opts_chunk$set(error=TRUE)
library(knitcitations)
library(nonparametricbayes)
library(reml)
```

```{r}
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

```

Rather than attempt a large list of nested loops, we do better to launch seperate scripts and write to a common data format.  We define a function that takes as arguments all things we wish to vary in the senstivity analysis, and returns a data.frame with the values of interest

```{r }
sensitivity <- function(model, parameters, nuisance, seed){

  if(model == "Myers")
    f <- Myers
  else if(model == "Allen")
    f <- RickerAllee

sigma_g <- nuisance[["sigma_g"]]
p <- parameters
sigma_m <- 0.0
z_g <- function() rlnorm(1, 0, sigma_g)
z_m <- function() 1
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)

                     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])

                     ## @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 mle-opt
                     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)

                     OPT = data.frame(GP = opt_gp$D, True = opt_true$D)
                     policies <- melt(data.frame(stock=x_grid, sapply(OPT, function(x) x_grid[x])), id="stock")
                     names(policies) <- c("stock", "method", "value")

                     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]]))
                     sims_data <- data.table(dat)
                     setnames(sims_data, c("L1", "L2"), c("method", "reps"))
                     # Legend in original ordering please, not alphabetical:
                     sims_data$method = factor(sims_data$method, ordered=TRUE, levels=names(OPT))

                     Profit <- sims_data[, sum(profit), by=c("reps", "method")]
                     tmp <- dcast(Profit, reps ~ method)
                     tmp <- tmp / tmp[,"True"]
                     tmp <- melt(tmp[2:dim(tmp)[2]])
                     tmp$value[tmp$variable == "GP"]
                     })

  yields_dat <- melt(yields)
  names(yields_dat) <- c("replicate", "simulation", "value")

  # Make ids factors, not integers
  yields_dat$replicate <- as.factor(yields_dat$replicate)
  yields_dat$simulation <- as.factor(yields_dat$simulation)

  # definitions of id codes.  In this case the id number is it's own definition.
  rep_ids <-  levels(yields_dat$replicate)
  names(rep_ids) <- rep_ids
  sim_ids <-  levels(yields_dat$simulation)
  names(sim_ids) <- sim_ids

  dat <- data.frame(model = model,
                    pars = as.list(parameters),
                    replicate = yields_dat$replicate,
                    sim = yields_dat$simulation,
                    value = yields_dat$value,
                    noise = sigma_g)
}
```

```{r}
model <- "Allen"
allen1.01 <- sensitivity(model,
                   parameters = parameters[[model]][[1]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
                   seed=c(1234, 2222, 3333))
```

```{r}
model <- "Allen"
allen2.01 <- sensitivity(model,
                   parameters = parameters[[model]][[2]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
                   seed=c(1234, 2222, 3333))
```

```{r}
model <- "Allen"
allen1.05 <- sensitivity(model,
                   parameters = parameters[[model]][[1]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
                   seed=c(1234, 2222, 3333))
```

```{r}
model <- "Allen"
allen2.05 <- sensitivity(model,
                   parameters = parameters[[model]][[2]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
                   seed=c(1234, 2222, 3333))
```

```{r}
model <- "Myers"
Myers1.01 <- sensitivity(model,
                   parameters = parameters[[model]][[1]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
                   seed=c(1234, 2222, 3333))
```

```{r}
model <- "Myers"
Myers2.01 <- sensitivity(model,
                   parameters = parameters[[model]][[2]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
                   seed=c(1234, 2222, 3333))
```

```{r}
model <- "Myers"
Myers1.05 <- sensitivity(model,
                   parameters = parameters[[model]][[1]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
                   seed=c(1234, 2222, 3333))
```

```{r}
model <- "Myers"
Myers2.05 <- sensitivity(model,
                   parameters = parameters[[model]][[2]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
                   seed=c(1234, 2222, 3333))
```

```{r}
save(list=ls(), file="sensitivity.rda")
```

```{r data}
allen_dat <- rbind(allen1.01, allen1.05,
             allen2.01, allen2.05)
m <- rbind(Myers1.01, Myers1.05,
             Myers2.01, Myers2.05)
myers_dat <- m[c(1:2,4,3,5:8)]
names(myers_dat) <- names(allen_dat)
model_dat <- rbind(allen_dat, myers_dat)
```

```{r}
dat <- model_dat
dat$pars.r <- factor(dat$pars.r, labels=c("A", "B", "C", "D"))
dat <- dat[c(1:2,5:6, 8, 7)]
dat$noise <- factor(dat$noise)
names(dat) <- c("model", "parameters", "replicate", "simulation", "noise", "value")
```

```{r}
ggplot(dat) + geom_histogram(aes(value)) + xlim(0,1.0) + theme_bw() + xlab("value as fraction of the optimal")

ggplot(dat) + geom_histogram(aes(value, fill=noise)) + xlim(0,1.0) + theme_bw() + xlab("value as fraction of the optimal")
```

```{r sensitivity}
ggplot(dat) + geom_histogram(aes(value, fill=noise)) + xlim(0,1.0) + theme_bw() + xlab("value as fraction of the optimal") + facet_wrap(~model, ncol=1)
```

```{r}
## Extract acutal parameter values corresponding to each parameter set
  p1 = levels(factor(model_dat$pars.r))
  p2 = levels(factor(model_dat$pars.K))
  p3 = levels(factor(model_dat$pars.C))

  A = c(r = p1[1], K = p2[1], theta = p3[1])
  B = c(r = p1[2], K = p2[2], theta = p3[2])
  C = c(r = p1[3], K = p2[3], C = p3[3])
  D = c(r = p1[4], K = p2[4], C = p3[4])
```

```{r eml, dependson='data'}
  rep_ids <-  levels(dat$replicate)
  sim_ids <-  levels(dat$sim)

  require(reml)
  col.defs = c(model = "name of the model used to simulate the data",
               parameters = "which set of model parameters was used",
               replicate = "replicate id numbers.  Each unique simulation of the model results in a policy function which is then replicated under each replicate id number to determine the realized value relative to the optimal value",
               simulation = "simulation id number. For each set of parameters for the model, multiple simulations of training data are generated and assigned unique simulation numbers",
               noise = "sigma_g value, multiplicative growth noise scale, as a log-normal log-standard deviation"
               value = "the value derived from this realization of the policy function estimated with the Gaussian process relative to the theoretical optimum value if the true model and parameter values were known",)

  unit.defs = list(model = c("Myers" = "a 3 parameter model of Beverton-Holt-like recruitment with the addition of an Allee effect (and fourth parameter for scaling the log-normal growth noise); based on a Myers et al 1995 (doi:10.1126/science.269.5227.1106)",
                               "Allen" = "a 3 parameter model of Ricker-like recruitment with the addition of an Allee effect (and fourth parameter for scaling the log-normal growth noise); based on Allen et al 2005 (doi:10.1080/10236190412331335373)"),
                  parameters =
  c("A" = paste("First set of parameter values drawn at random for Myers model. r, K, theta respectively are:", A),
    "B" = paste("Second set of parameter values drawn at random for Myers model. r, K, theta respectively are:", B),
    "C" = paste("First set of parameter values drawn at random for Myers model. r, K, C respectively are:", C),
    "D" = paste("Second set of parameter values drawn at random for Myers model. r, K, C respectively are:", D)),
                  replicate = rep_ids,
                  simulation = sim_ids,
                  noise = "density/density (pure)",
                  value = "Realized Dollars/Potential Dollars")

emldat =  data.set(dat, col.defs=col.defs, unit.defs=unit.defs)
eml_write(emldat, creator="Carl Boettiger <cboettig@gmail.com>")
```

```{r, eval=FALSE}
require(reshape2)
require(snowfall)
sfInit(parallel=TRUE, cpu=8)
dat <- sfLapply(nuisance_values[["sigma_g"]], function(sigma){
  dat <- lapply(parameters[["Allen"]], function(p){
         sensitivity("Allen",
                   parameters = p,
                   nuisance = c(sigma_g = sigma),
                   seed=c(1234, 2222, 3333))
              })
  dat <- melt(dat, id=names(dat[[1]]))
})
dat <- melt(dat, id=names(dat[[1]]))
```
```{r eval=FALSE}
require(reshape2)
myers_dat <- sfLapply(nuisance_values[["sigma_g"]], function(sigma){
  dat <- lapply(parameters[["Myers"]], function(p){
         sensitivity("Myers",
                   parameters = p,
                   nuisance = c(sigma_g = sigma),
                   seed=c(1234, 2222, 3333))
              })
  dat <- melt(dat, id=names(dat[[1]]))
})
myers_dat <- melt(myers_dat, id=names(dat[[1]]))

```


text:
require(knitr)
knit("sensitivity.Rmd")


text:
```{r}
## 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)
```

```{r}
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
}
```

```{r}
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]]))
}
```

```{r}
model <- "Allen"
allen1.01 <- sensitivity(model,
                   parameters = parameters[[model]][[1]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
                   seed=c(1234, 2222, 3333))
```
```{r}
model <- "Allen"
allen2.01 <- sensitivity(model,
                   parameters = parameters[[model]][[2]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
                   seed=c(1234, 2222, 3333))
```
```{r}
model <- "Allen"
allen1.05 <- sensitivity(model,
                   parameters = parameters[[model]][[1]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
                   seed=c(1234, 2222, 3333))
```
```{r}
model <- "Allen"
allen2.05 <- sensitivity(model,
                   parameters = parameters[[model]][[2]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
                   seed=c(1234, 2222, 3333))

```
```{r}
model <- "Myers"
Myers1.01 <- sensitivity(model,
                   parameters = parameters[[model]][[1]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
                   seed=c(1234, 2222, 3333))
```
```{r}
model <- "Myers"
Myers2.01 <- sensitivity(model,
                   parameters = parameters[[model]][[2]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[1]),
                   seed=c(1234, 2222, 3333))
```
```{r}
model <- "Myers"
Myers1.05 <- sensitivity(model,
                   parameters = parameters[[model]][[1]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
                   seed=c(1234, 2222, 3333))

```
```{r}
model <- "Myers"
Myers2.05 <- sensitivity(model,
                   parameters = parameters[[model]][[2]],
                   nuisance = c(sigma_g = nuisance_values$sigma_g[2]),
                   seed=c(1234, 2222, 3333))
```

```{r}
## 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)
```


text:
# Learning ADMB

Plotting and knitr options, (can generally be ignored)

```{r plotting-options, cache=FALSE, message=FALSE, warning=FALSE, include=FALSE}
source("~/.knitr_defaults.R")
opts_chunk$set(warning=TRUE, message=TRUE, cache=FALSE)
opts_knit$set(upload.fun = socialR::flickr.url)
```

### Model and parameters

```{r stateeq}
f <- function(x,h,p)  x * exp(p[1] * (1 - x / p[2]) * (x - p[3]) / p[2] )
p <- c(1, 10, 5)
K <- 10  # approx, a li'l' less
Xo <- 6 # approx, a li'l' less
```

Various parameters defining noise dynamics, grid, and policy costs.

```{r sdp-pars, dependson="stateeq"}
sigma_g <- 0.1
z_g <- function() rlnorm(1,0, sigma_g)
x_grid <- seq(0, 1.5 * K, length=50)
Tobs <- 40
set.seed(123)
```

### Sample Data

```{r obs, dependson="sdp-pars"}
x <- numeric(Tobs)
x[1] <- Xo
for(t in 1:(Tobs-1))
  x[t+1] = z_g() * f(x[t], h=0, p=p)
qplot(1:Tobs, x)
```

## Maximum Likelihood "by hand"

```{r mle, dependson="obs"}
STABLIZE = 1e-10
n = length(x)
mloglik <- function(pars){
  r = pars[1]; k = pars[2]; c = pars[3]; s = pars[4];
  mu = (x+STABLIZE) * exp( r * (1 - x / (k+STABLIZE)) * (x - c) / (k + STABLIZE));
  mu = pmin(1e100, mu) # avoid infinite values
  f = 0.5 * n * log(2 * pi) + n * log(s + STABLIZE) + 0.5 * sum(x - mu + STABLIZE)^2/ (s + STABLIZE)^2;

  f
  }
```

Starting from the true values we mostly just shrink the noise parameter:

```{r}
init <- c(p, sigma_g)
mloglik(init) #true minus loglik
o <- optim(init, mloglik, method="L", lower=1e-5, upper=1e5)
o$value
o$par
```

While starting from arbitrary values we still find the optim.

```{r}
init <- c(1,1,1,1)
init <- c(p, sigma_g)
mloglik(init) #true minus loglik
o <- optim(init, mloglik, method="L", lower=1e-5, upper=1e5)
o$value
o$par
```

Okay, now lets try admb.  We use R2admb which is just a convenient way to write our data and parameters into an admb file.

```{r admb-lib}
# install_github("R2admb", "bbolker", subdir="R2admb") # dev version
library(R2admb)
```

## ADMB definition

We still need to define the model using ADMB notation in the procedure section.  This is mostly like R or C++, with the exception of special functions like `square` in place of `^2`, `norm2` for the sum of squares, and `elem_prod` istead of `*` for the element-wise product of two arrays. The constant `pi` is given as `M_PI`, as typical of C/C++ libraries.  Where these other functions are defined I'm not sure, but some useful guides to [ADMB vector/matrix operations](http://fish.washington.edu/research/MPAM/resources/ADMB_Minte-Vera.pdf) or an (undefined) list of [keywords](http://www.admb-project.org/developers/contribute-documentation/functions/keywords.txt/view)...

The equivalent model

```{r model}
model <-
paste("
PARAMETER_SECTION
  vector mu(1,n) // per capita mort prob

PROCEDURE_SECTION
  mu = log(x) + r * elem_prod((1 - x / k), (x - c) / k);
  f = 0.5 * n * log(2 * M_PI) + n * log(s) + 0.5 * norm2(x - exp(mu)) / square(s);
")
writeLines(model, "model.tpl")
```

Without explicit handling of the overflow errors, ADMB does not give us reliable estimates with arbitrary starting conditions

```{r run-admb1}
setup_admb("/var/admb")

df <- data.frame(x=x)
params <- list(r = 1, k = 1, c = 1, s = 1) ## starting parameters
bounds <- list(r = c(1e-10, 1e3), k=c(1e-10, 1e3), c=c(1e-10, 1e3), s = c(1e-5,1e3)) ## bounds
dat <- c(list(n = nrow(df)), df)
m1 <- do_admb("model",
              data = dat,
              params = params,
              bounds = bounds,
              run.opts = run.control(checkparam="write",
                                     checkdata="write", clean=FALSE))
m1
```

But does fine with good starting values.  Hmm.. thought that was supposed to be the other way around...

```{r run-admb}
params <- list(r = 1, k = 10, c = 5, s = .1) ## starting parameters

m1 <- do_admb("model",
              data = dat,
              params = params,
              bounds = bounds,
              run.opts = run.control(checkparam="write",
                                     checkdata="write"))
m1
```

Which finds a better optim (though substantailly overfit in reality)

Hans suggests adding an error term in the function definitions rather than in limiting the bounds or log transforming the variables:

> The most common plase where
 this goes wrong is 1/0, log(0), sqrt(0), pow(0,1) etc.
> Your suggestion is OK, but usually I prefer to put
 in log(1e-10+my_expression), sqrt(1e-10+my_expression), pow(1e-10+my_expression,1)


text:
```{r set-options, echo = FALSE, cache = FALSE, external = TRUE, include = FALSE}
opts_chunk$set(external = TRUE, cache = FALSE, cache.path = "myers-cache/", warning=FALSE)
read_chunk('gaussian-process-control.R')
library(knitcitations)
```

```{r libraries, include=FALSE}
```
```{r graphing-options, include=FALSE}
```

```{r Myer}
```

We use the model of `r citet("10.1126/science.269.5227.1106")`.

```{r sdp-pars}
```

With parameters `r p`.

```{r x0-obs}
x_0_observed <- allee + x_grid[5]
xT <- 0
set.seed(1)
```

```{r sim-obs}
```

We simulate data under this model, starting from a size of `r x_0_observed`.

```{r lag-data}
```

We consider the observations as ordered pairs of observations of current stock size $x_t$ and observed stock in the following year, $x_{t+1}$.  We add the pseudo-observation of $0,0$.  Alternatively we could condition strictly on solutions passing through the origin, though in practice the weaker assumption is often sufficient.

```{r par-est}
```

```{r gp-priors}
```

```{r gp-fit, dependson="lag-data"}
```

We fit a Gaussian process with

```{r gp-data, dependson="gp-fit"}
```

```{r gp-plot, dependson="gp-data"}
```

```{r gp-posteriors}
```

## The transition matrix of the inferred process

```{r persistence-test}
```

```{r gp-F-sim, dependson=c("persistence-test", "gp-data")}
```

```{r par-F-sim, dependson="persistence-test"}
```

```{r F-sim-plot, dependson=c("gp-F-sim", "par-F-sim")}
```

```{r gp-opt, dependson="gp-data"}
```

```{r true-opt}
```

```{r est-opt}
```

```{r policy_plot}
```

```{r othernoise}
```

```{r stationary_policy_only}
```

```{r simulate}
```

```{r tidy}
```

```{r sim-fish}
```

```{r sim-harvest}
```

```{r costs}
```

```{r echo=FALSE, results="asis"}
bibliography("html")
```


text:
DATA_SECTION

  init_int n
  init_vector x(1,40)

PARAMETER_SECTION

  objective_function_value f
  init_bounded_number r(0,100)
  init_bounded_number k(0,100)
  init_bounded_number c(0,100)
  init_bounded_number s(0,100)
  vector mu(1,n) // per capita mort prob
PROCEDURE_SECTION
  mu = x + r * elem_prod((1 - x / k), (x - c) / k);
  f = 0.5 * n * log(2 * M_PI) + n * log(s) + 0.5 * norm2(x - mu) / square(s);


text:
# Semi-analytic posteriors

### Generating Model and parameters

Ricker model, parameterized as

$$X_{t+1} = X_t r e^{-\beta X_t + \sigma Z_t}$$

for random unit normal $Z_t$

```r
f <- function(x,h,p)  x * p[1] * exp(-x * p[2])
p <- c(exp(1), 1/10)
K <- 10  # approx, a li'l' less
Xo <- 1 # approx, a li'l' less
```

```r
sigma_g <- 0.1
z_g <- function() rlnorm(1,0, sigma_g)
x_grid <- seq(0, 1.5 * K, length=50)
N <- 40
set.seed(123)
```

### Sample Data

```r
x <- numeric(N)
x[1] <- Xo
for(t in 1:(N-1))
  x[t+1] = z_g() * f(x[t], h=0, p=p)
qplot(1:N, x)
```

![plot of chunk obs](http://farm4.staticflickr.com/3747/8965288320_9c74f855e3_o.png)

## Compute the posterior after marginalizing over $r$ and $\sigma$ parameters:

$$P(\beta | X) $$

```r
Mt <- function(t, beta) log(x[t+1]) - log(x[t]) + beta * x[t]
beta_grid = seq(1e-5, 2, by=1e-3)

P_B.X <- sapply(beta_grid, function(beta){
  Mt_vec = sapply(1:(N-1), Mt, beta)
  sum_of_squares <- sum(Mt_vec^2)
  square_of_sums <- sum(Mt_vec)^2
  0.5 ^ (N/2-1) * (sum_of_squares - square_of_sums/(N-1)) ^ (N/2-1) / gamma(N/2-1)
  })

qplot(beta_grid, -log(P_B.X))
```

![plot of chunk unnamed-chunk-1](http://farm3.staticflickr.com/2852/8965288614_24b1714461_o.png)

```r

beta_grid[which.min(P_B.X)]
```

```
[1] 0.09801
```

Estimating the Allen model:

$$X_{t+1} = X_t \exp\left(r \left(1 - \frac{X_t}{K} \right\left(\frac{X - C}{K}\right)\right) $$

First re-parameterize so we can isolate an additive parameter, using standard quadratic form,

$$X_{t+1} = X_t \exp(c + b X_t + a X_t^2) $$

Where $c = \tfrac{-rC}{K}$, $b = \tfrac{r}{K}\left(\tfrac{C}{K} + 1\right)$, and $a = \tfrac{r}{K^2}$.  Then after integrating out $c$ and $\sigma$ we have

```r
Mt <- function(t, a, b)  log(x[t+1]) - log(x[t]) - b * x[t] + a * x[t] ^ 2
```

Choosing a sensible grid is a bit more tricky following the transformation, particularly $b$ can be negative or positive, and much larger in magnitude than any of the original parameters.  $a$ on the other hand, can be reasonably restricted:

```r
b_grid = seq(0, .1, length=100)
a_grid = seq(0, .1, length=100)
```

We define the probability as a function of the remaining two parameters,

```r
prob <- function(a, b){
  Mt_vec = sapply(1:(N-1), Mt, a, b)
  sum_of_squares <- sum(Mt_vec^2)
  square_of_sums <- sum(Mt_vec)^2
0.5 ^ (N/2-1) * (sum_of_squares - square_of_sums/(N-1)) ^ (N/2-1) / gamma(N/2-1)
}
```

and loop over the grid on each

```r
P_a_b.X <- sapply(a_grid, function(a)
                sapply(b_grid, function(b) prob(a, b)))
```

We summarize with a contour plot

```r
  df = melt(P_a_b.X)
names(df) = c("a", "b", "lik")
ggplot(df, aes(a_grid[a], b_grid[b], z=-log(lik))) +
  stat_contour(aes(color=..level..), binwidth=1) +
  coord_cartesian(ylim=c(0, .025), xlim=c(0,.05))
```

![plot of chunk unnamed-chunk-6](http://farm9.staticflickr.com/8125/8964095361_4bf8fa162d_o.png)

Estimating the Myers model on this data:

$$X_{t+1} = Z_t \frac{r X_t^{\theta}}{1 + X_t^{\theta} / K}$$

With $Z_t$ lognormal, unit mean, std $\sigma$.

```r
Mt <- function(t, theta, K) log(x[t+1]) - theta * log(x[t]) + log(1 + x[t] ^ theta / K)
theta_grid = seq(1e-5, 5, length=100)
K_grid = seq(1e-5, 30, length=100)

prob <- function(theta, K){
  Mt_vec = sapply(1:(N-1), Mt, theta, K)
  sum_of_squares <- sum(Mt_vec^2)
  square_of_sums <- sum(Mt_vec)^2
0.5 ^ (N/2-1) * (sum_of_squares - square_of_sums/(N-1)) ^ (N/2-1) / gamma(N/2-1)
}

P_theta_K.X <- sapply(theta_grid, function(theta)
                sapply(K_grid, function(k) prob(theta, k)))

require(reshape2)
df = melt(P_theta_K.X)
names(df) = c("theta", "K", "lik")
ggplot(df, aes(theta_grid[theta], K_grid[K], z=-log(lik))) + stat_contour(aes(color=..level..), binwidth=3)
```

![plot of chunk unnamed-chunk-8](http://farm4.staticflickr.com/3690/8964097371_85d0c17d08_o.png)


text:
Fixed priors on hyperparameters, fixed model type.

```r
#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)
tau2.p <- c(5,1)
d.p = c(10, 1/0.1, 10, 1/0.1)
nug.p = c(10, 1/0.1, 10, 1/0.1) # gamma mean
s2_prior <- function(x) dinvgamma(x, s2.p[1], s2.p[2])
tau2_prior <- function(x) dinvgamma(x, tau2.p[1], tau2.p[2])
d_prior <- function(x) dgamma(x, d.p[1], scale = d.p[2]) + dgamma(x, d.p[3], scale = d.p[4])
nug_prior <- function(x) dgamma(x, nug.p[1], scale = nug.p[2]) + dgamma(x, nug.p[3], scale = nug.p[4])
beta0_prior <- function(x, tau) dnorm(x, 0, tau)
beta = c(0)
priors <- list(s2 = s2_prior, tau2 = tau2_prior, beta0 = dnorm, nug = nug_prior, d = d_prior, ldetK = function(x) 0)
```

```r
profit = function(x,h) pmin(x, h)
delta <- 0.01
OptTime = 20  # stationarity with unstable models is tricky thing
reward = 0
xT <- 0
```

```r
f <- RickerAllee
# c(5, 10, 5) is 2-cycle, c(5.5, 10, 5) is 6 cycle, 5.3 is about 4
p <- c(2, 10, 5)
K <- 10
allee <- 5
```

```r
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=101)
h_grid <- x_grid
```

With parameters 2, 10, 5.

```r
seed_i <- 1
  Xo <- K # observations start from
  x0 <- Xo # simulation under policy starts from
```

```r
  obs <- sim_obs(Xo, z_g, f, p, Tobs=35, nz=15,
                 harvest = sort(rep(seq(0, .5, length=7), 5)), seed = seed_i)
```

```r
  alt <- par_est(obs,  init = c(r=p[1], K=mean(obs$x), s=sigma_g))
  est <- par_est_allee(obs, f, p,  init = c(2, mean(obs$x), 2, s = sigma_g))
```

Which estimates a Ricker model with $r =$ 2, $K =$ 6.77, and the Allen allee model with $r =$ 2, $K =$ 6.77 and $C =$ 2.

```r
  gp <- bgp(X=obs$x, XX=x_grid, Z=obs$y, verb=0,
          meanfn="constant", bprior="b0", BTE=c(10,1600,2),
          m0r1=FALSE, corr="exp", trace=TRUE,
          beta = beta, s2.p = s2.p, d.p = d.p, nug.p = nug.p, tau2.p = tau2.p,
          s2.lam = "fixed", d.lam = "fixed", nug.lam = "fixed", tau2.lam = "fixed")
  gp_plot(gp, f, p, est$f, est$p, alt$f, alt$p, x_grid, obs, seed_i)
```

![plot of chunk unnamed-chunk-4](http://farm9.staticflickr.com/8401/8693278143_1b8d80d884_o.png)

```r
  posteriors_plot(gp, priors) # needs trace=TRUE!
```

```
stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
this.
```

```
stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
this.
```

```
stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
this.
```

```
stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
this.
```

```
stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
this.
```

```
stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust
this.
```

![plot of chunk unnamed-chunk-5](http://farm9.staticflickr.com/8254/8693278361_6a8c45c339_o.png)

```r
  OPT <- optimal_policy(gp, f, est$f, alt$f,
                        p, est$p, alt$p,
                        x_grid, h_grid, sigma_g,
                        sigma_g, sigma_g, # est$sigma_g, alt$sigma_g, but those ests are poor
                        delta, xT, profit, reward, OptTime)
  plot_policies(x_grid, OPT$gp_D, OPT$est_D, OPT$true_D, OPT$alt_D)
```

![plot of chunk unnamed-chunk-6](http://farm9.staticflickr.com/8260/8693286775_966faae4ca_o.png)

```r
dt <- simulate_opt(OPT, f, p, x_grid, h_grid, x0, z_g, profit, OptTime=OptTime)
sim_plots(dt, seed=seed_i)
```

```
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
```

![plot of chunk unnamed-chunk-7](http://farm9.staticflickr.com/8118/8693287195_e1ff8acbc6_o.png)

```r
profits_stats(dt)
```

```
       method     V1     sd
1:         GP 19.337 1.6314
2: Parametric  6.446 0.3384
3:       True 20.701 1.8012
4: Structural  7.650 0.0000
```