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}
require(pdgControl)
require(nonparametricbayes)
require(reshape2)
require(ggplot2)
require(data.table)
require(tgp)
library(kernlab)
require(MCMCpack)
require(plyr)

```

```{r graphing-options, include=FALSE}
opts_knit$set(upload.fun = socialR::notebook.url)
#opts_chunk$set(dev = 'Cairo_pdf', dev.args=list(""))
opts_chunk$set(dev="png", dev.args=list(bg="transparent"))
opts_chunk$set(comment=NA, tidy=FALSE)
theme_set(theme_bw(base_size=16))
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")
```

Fixed priors on hyperparameters, fixed model type.

```{r gp-priors}
s2.p <- c(50,50)  # inverse gamma mean is
tau2.p <- c(20,1)
d.p = c(10, 1/0.01, 10, 1/0.01)
nug.p = c(10, 1/0.01, 10, 1/0.01) # 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 sdp-pars-fixed}
profit = function(x,h) pmin(x, h)
delta <- 0.01
OptTime = 20  # stationarity with unstable models is tricky thing
reward = 0
xT <- 0
z_g = function() rlnorm(1, 0, sigma_g)
z_m = function() 1+(2*runif(1, 0,  1)-1) * sigma_m
```

```{r RickerAllee-exp}
f <- RickerAllee
p <- c(1.1, 10, 5)
K <- 10
allee <- 5
```

```{r sdp-pars-explore}
sigma_g <- 0.05
sigma_m <- 0.02
x_grid <- seq(0, 1.5 * K, length=101)
h_grid <- x_grid
```

With parameters `r p`.

```{r}

```

```{r parallel, message=FALSE, warning=FALSE, results='hide'}
require(snowfall)
sfInit(par=TRUE, cpu=8)
sfExportAll()
sfLibrary(pdgControl)
sfLibrary(nonparametricbayes)
sfLibrary(reshape2)
sfLibrary(ggplot2)
sfLibrary(data.table)
sfLibrary(tgp)
sfLibrary(kernlab)
sfLibrary(MCMCpack)
sfLibrary(plyr)
```

```{r }
seed <- round(runif(32) * 1e6)
seed[1] <- 1
seed
yields <- lapply(seed,
function(seed_i){
  Xo <- allee + x_grid[10] # observations start from
  x0 <- Xo # simulation under policy starts from
  obs <- sim_obs(Xo, z_g, f, p, Tobs=40, seed = seed_i)
  est <- par_est(obs)
  gp <- bgp(X=obs$x, XX=x_grid, Z=obs$y, verb=0,
          meanfn="constant", bprior="b0", BTE=c(2000,16000,2),
          m0r1=FALSE, corr="exp", trace=FALSE,
          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_alt, est$p_alt, x_grid, obs, seed_i)
#  posteriors_plot(gp, priors) # needs trace=TRUE!
  OPT <- optimal_policy(gp, f, est$f_alt, p, est$p_alt, x_grid, h_grid, sigma_g, est$sigma_g_alt, delta, xT, profit, reward, OptTime)
  dt <- simulate_opt(OPT, f, p, x_grid, h_grid, x0, z_g, profit)
  sim_plots(dt, seed=seed_i)
  profits_stats(dt)
})
```

```{r}
yields <- melt(yields, id=c("method", "V1", "sd"))
ggplot(yields) + geom_histogram(aes(V1)) + facet_wrap(~method, scale='free')
yields
````

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


text:
```{r plotting-options}
opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE, cache=FALSE, comment=NA, fig.width=6, fig.height=4)
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))
cbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
               "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
```

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

Required R libraries (for multivariate normal, also for plotting):

```{r results="hide"}
require(MASS)
require(reshape2)
```

Set a seed for repeatable plots

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

Define the points at which we want to compute the function values (x values of the prediction points or test points), and the scale parameter for the covariance function $\ell=1$

```{r}
x_predict <- seq(-5,5,len=50)
l <- 1
```

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

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

```{r}
COV <- cov(x_predict, x_predict)
```
Generate a number of functions from the process

```{r}
values <- mvrnorm(3, rep(0, length=length(x_predict)), COV)
```

Reshape the data into long (tidy) form, listing x value, y value, and sample number
```{r}
dat <- data.frame(x=x_predict, t(values))
dat <- melt(dat, id="x")
head(dat)
```

Plot the result

```{r fig2a}
fig2a <- ggplot(dat,aes(x=x,y=value)) +
  geom_rect(xmin=-Inf, xmax=Inf, ymin=-2, ymax=2, fill="grey80") +
  geom_line(aes(group=variable)) +
  scale_y_continuous(lim=c(-3,3), name="output, f(x)") +
  xlab("input, x")
fig2a
```

### Posterior distribution given the data

Imagine we have some data,

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

In general we aren't interested in drawing from the prior, but want to include information from the data as well.  We want the joint distribution of the observed values and the prior is:

$$\begin{pmatrix} y_{\textrm{obs}} \\ y_{\textrm{pred}} \end{pmatrix} \sim \mathcal{N}\left( \mathbf{0}, \begin{bmatrix} cov(X_o,X_o) & cov(X_o, X_p) \\ cov(X_p,X_o) & cov(X_p, X_p) \end{bmatrix} \right)$$

## No observation noise

Assuming the data are known without error and conditioning on the data, and given $x \sim \mathcal{N}(0, cov(X_o, X_o))$, then the conditional probability of observing our data is easily solved by exploiting the nice properties of Gaussians,

$$x|y \sim \mathcal{N}(E,C)$$
$$E = cov(X_p, X_o) cov(X_o,X_o)^{-1} y$$
$$C= cov(X_p, X_p) - cov(X_p, X_o) cov(X_o,X_o)^{-1} cov(X_o, X_p)$$

(We use `solve(M)` which with no second argument will simply invert the matrix `M`, but should use the Cholsky decomposition instead for better numerical stability)

```{r}
cov_xx_inv <- solve(cov(obs$x, 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)
```

Now lets take 3 random samples from the posterior distribution,
```{r}
values <- mvrnorm(3, Ef, Cf)
```

and plot our solution (mean, 2 standard deviations, and the ranom samples.)

```{r fig2b}
dat <- data.frame(x=x_predict, t(values))
dat <- melt(dat, id="x")

fig2b <- ggplot(dat,aes(x=x,y=value)) +
  geom_ribbon(data=NULL,
              aes(x=x_predict, y=Ef, ymin=(Ef-2*sqrt(diag(Cf))), ymax=(Ef+2*sqrt(diag(Cf)))),
              fill="grey80") +
  geom_line(aes(color=variable)) + #REPLICATES
  geom_line(data=NULL,aes(x=x_predict,y=Ef), size=1) + #MEAN
  geom_point(data=obs,aes(x=x,y=y)) +  #OBSERVED DATA
  scale_y_continuous(lim=c(-3,3), name="output, f(x)") +
  xlab("input, x")
fig2b
```

```{r fig2_data}
ggplot(data=obs,aes(x=x,y=y)) +
  geom_point()   #OBSERVED DATA
```

Additive noise
--------------

In general the model may have process error, and rather than observe the deterministic map $f(x)$ we only observe $y = f(x) + \varepsilon$.  Let us assume for the moment that $\varepsilon$ are independent, normally distributed random variables with variance $\sigma_n^2$.

```{r}
sigma.n <- 0.3
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)
```

Now lets take 3 random samples from the posterior distribution,
```{r}
values <- mvrnorm(3, Ef, Cf)
```

and plot

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

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

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

repo_name: cboettig/nonparametric-bayes

repo_language: TeX

mime_type: text/plain


text:
Bayesian regression in Gaussian processes
---------------------------------------

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

![plot of chunk unnamed-chunk-6](http://farm8.staticflickr.com/7285/8739390350_c4c06db5ab_o.png) ![plot of chunk unnamed-chunk-6](http://farm8.staticflickr.com/7285/8739390438_50350c1a38_o.png)

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

![plot of chunk unnamed-chunk-8](http://farm8.staticflickr.com/7282/8738272607_d913f04a10_o.png)


text:
## @knitr libraries
require(pdgControl)
require(nonparametricbayes)
require(reshape2)
require(ggplot2)
require(data.table)
require(tgp)
library(kernlab)
require(MCMCpack)
require(plyr)

## @knitr graphing-options
opts_knit$set(upload.fun = function(file) socialR::notebook.url(file, cp=TRUE, sync=FALSE))
#opts_chunk$set(dev = 'Cairo_pdf', dev.args=list(""))
opts_chunk$set(dev="png", dev.args=list(bg="transparent"))
opts_chunk$set(comment=NA, tidy=FALSE)
theme_set(theme_bw(base_size=16))
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")
#sha <- gsub("^commit ", "", system("git log -n 1", intern=TRUE)[1])
#short_sha <- gsub("(^.{10}).*", "\\1", sha)
#date <- format(Sys.time(), "%Y-%m-%d-%H-%M-%S")
#opts_chunk$set(fig.path = paste("figure/", date, "-", short_sha, "-", sep=""))

## @knitr BevHolt
f <- BevHolt
p <- c(1.5,.05)
K <- (p[1]-1)/p[2]

## @knitr Myer
f <- Myer_harvest
pars <- c(1, 2, 4.5)
p <- pars # shorthand
K <- p[1] * p[3] / 2 + sqrt( (p[1] * p[3]) ^ 2 - 4 * p[3] ) / 2
allee <- p[1] * p[3] / 2 - sqrt( (p[1] * p[3]) ^ 2 - 4 * p[3] ) / 2 # allee threshold
e_star <- (p[1] * sqrt(p[3]) - 2) / 2 ## Bifurcation point

## @knitr Ricker
f <- Ricker
p <- c(2.5, 10)
K <- 12.5

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

## @knitr plot-may
birth <- function(x) p["r"] * (1-  x / p["k"])
death <- function(x) p["a"] * x ^ (p["Q"] - 1) / (x ^ p["Q"] + p["H"])
df <- data.frame(x = x_grid, b = sapply(x_grid, birth), d = sapply(x_grid, death))
ggplot(df) + geom_line(aes(x, b), col = "blue") + geom_line(aes(x,d), col = "red")

## @knitr RickerAllee
f <- RickerAllee
p <- c(r = 3, K = 10, C=4)
K <- p[2]
allee <- p[3]

## @knitr sdp-pars
sigma_g <- 0.05
sigma_m <- 0.2
z_g <- function() 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 = 50 # should be long enough that profit is greater sustainably harvesting
reward = profit(x_grid[length(x_grid)], x_grid[length(x_grid)]) + 1 / (1 - delta) ^ OptTime
xT <- 0

## @knitr inits
# x_0_observed is starting condition for simulation of the observed data.
# It should be in preferred state for bistable model,
# above Allee threshold for Allee model,
# and near zero for BH or Ricker models
xT <- allee
x_0_observed <- K

## @knitr sim-obs
Tobs <- 30
x <- numeric(Tobs)
x[1] <- x_0_observed
for(t in 1:(Tobs-1))
  x[t+1] = z_g() * f(x[t], h=0, p=p)
plot(x)

## @knitr lag-data
obs <- data.frame(x=c(0,x[1:(Tobs-1)]),y=c(0,x[2:Tobs]))

## @knitr sim-with-harvest
nz <- 10 # weight on 0,0 pt
Tobs <- 40
harvest <- sort(rep(seq(1, 2, length=5), 8))
x <- numeric(Tobs)
x[1] <- x_0_observed
for(t in 1:(Tobs-1))
  x[t+1] = z_g() * f(x[t], h=2, p=p)
plot(x)

obs <- data.frame(x=c(rep(0,nz), pmax(rep(0,Tobs-1),x[1:(Tobs-1)]-harvest[1:Tobs-1])),y=c(rep(0,nz),x[2:Tobs]))
plot(obs$x, obs$y)

## @knitr MLE-RickerAllee
estf <- function(p){
  mu <- f(obs$x,0,p)
  -sum(dlnorm(obs$y, log(mu), p["s"]), log=TRUE)
}
par = c(r = p[1] - 1,
        C = p[2] - 1,
        K = p[3] + 2,
        s = sigma_g + abs(rnorm(1,0,.1)))
o <- optim(par, estf, method="L", lower=c(1e-3,1e-3,1e-3, 1e-3))
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])
p_alt
sigma_g_alt

## @knitr par-est
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']

## @knitr par-est-BH
estf <- function(p){
  mu <- obs$x * p["A"] / (1+ obs$x * p["B"])
  -sum(dlnorm(obs$y, log(mu), p["s"]), log=TRUE)
}
o <- optim(par = c(A=2, B=1/mean(x),s=1), estf, method="L", lower=c(1e-3,1e-3,1e-3))
f_alt <- BevHolt
p_alt <- c(A=as.numeric(o$par['A']), B=as.numeric(o$par['B']))
sigma_g_alt <- o$par['s']

## @knitr gp-priors
s2.p <- c(50,50) #inv gamma has mean b / (a - 1) (assuming a>1) and variance b ^ 2 / ((a - 2) * (a - 1) ^ 2) (assuming a>2)
tau2.p <- c(50,50)
d.p = c(10, 1/0.01, 10, 1/0.01)  ## sum of gammas, shape a & rate b. gamma has mean a / b and variance a / b ^ 2
nug.p = c(10, 1/0.01, 10, 1/0.01)
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 = tau2.p[2] / (tau2.p[1] - 1)) dnorm(x, 0, tau)
beta = c(0)

## @knitr gp-fit
gp <- bgp(X=obs$x, XX=x_grid, Z=obs$y, verb=0,
          meanfn="constant", bprior="b0", BTE=c(2000,16000,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")

## @knitr gp-data
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),
                 ymin2 = gp$ZZ.mean - 1.96 * sqrt(gp$ZZ.vark),
                 ymax2 = gp$ZZ.mean + 1.96 * sqrt(gp$ZZ.vark) )

# ZZ.ks2, ZZ.vark, ZZ.s2

## @knitr gp-plot
true <- sapply(x_grid, f, 0, p)
est <- sapply(x_grid, f_alt, 0, p_alt)
models <- data.frame(x=x_grid, GP=tgp_dat$y, Parametric=est, True=true)
models <- melt(models, id="x")
names(models) <- c("x", "method", "value")
# plot
ggplot(tgp_dat)  +
  geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") +
  geom_ribbon(aes(x,y,ymin=ymin2,ymax=ymax2), fill="gray60") +
  geom_line(data=models, aes(x, value, col=method), lwd=2, alpha=0.8) +
  geom_point(data=obs, aes(x,y), alpha=0.8) +
  xlab(expression(X[t])) + ylab(expression(X[t+1])) +
  scale_colour_manual(values=cbPalette)

## @knitr gp-posteriors
hyperparameters <- c("index", "s2", "tau2", "beta0", "nug", "d", "ldetK")
posteriors <- melt(gp$trace$XX[[1]][,hyperparameters], id="index")
priors <- list(s2 = s2_prior, tau2 = tau2_prior, beta0 = beta0_prior, nug = nug_prior, d = d_prior, ldetK = function(x) 0)
prior_curves <- ddply(posteriors, "variable", function(dd){
  grid <- seq(min(dd$value), max(dd$value), length = 100)
  data.frame(value = grid, density = priors[[dd$variable[1]]](grid))
})
ggplot(posteriors) +
  #geom_density(aes(value), lwd=2) +
  geom_histogram(aes(x=value, y=..density..), lwd=2) +
  geom_line(data=prior_curves, aes(x=value, y=density), col="red", lwd=2) +
  facet_wrap(~ variable, scale="free")
ggplot(prior_curves) +
  geom_line(aes(x=value, y=density), col="red", lwd=2) +
  facet_wrap(~ variable, scale="free")
#ggplot(subset(posteriors, variable=="nug")) + geom_histogram(aes(x=value, y = ..density..), lwd=2) + stat_function(fun = nug_prior, col="red", lwd=2)
#ggplot(subset(posteriors, variable=="s2")) + geom_histogram(aes(x=value, y = ..density..), lwd=2) + stat_function(fun = s2_prior, col="red", lwd=2)

## @knitr persistence-test
X <- numeric(length(x_grid))
X[38] = 1
h <- 0

## @knitr gp-F-sim
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")

## @knitr par-F-sim
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")

## @knitr F-sim-plot
transition <- melt(data.frame(x = x_grid, gp = xt1[1,], parametric = yt1[1,]), id="x")
ggplot(transition) + geom_point(aes(x,value, col=variable))

## @knitr est-F-sim
F_est <- par_F(h, f_alt, p_alt, x_grid, sigma_g)
zt1 <- X %*% F_est
zt10 <- zt1
for(s in 1:OptTime)
  zt10 <- zt10 %*% F_est
qplot(x_grid, zt10[1,]) + geom_point(aes(y=zt1[1,]), col="grey")

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

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

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

## @knitr policy_plot
policies <- melt(data.frame(stock=x_grid,
                            GP = x_grid[opt_gp$D[,1]],
                            Parametric = x_grid[opt_estimated$D[,1]],
                            True = x_grid[opt_true$D[,1]]),
                  id="stock")
names(policies) <- c("stock", "method", "value")
policy_plot <- ggplot(policies, aes(stock, stock - value, color=method)) +
  geom_line(lwd=2, alpha=0.8) +
  xlab("stock size") + ylab("escapement")  +
  scale_colour_manual(values=cbPalette)
policy_plot

## @knitr othernoise
z_g = function() rlnorm(1, 0, sigma_g)
z_m = function() 1+(2*runif(1, 0,  1)-1) * sigma_m

## @knitr stationary_policy_only
m <- sapply(1:OptTime, function(i) opt_gp$D[,1])
opt_gp$D <- m
mm <- sapply(1:OptTime, function(i) opt_true$D[,1])
opt_true$D <- mm
mmm <- sapply(1:OptTime, function(i) opt_estimated$D[,1])
opt_estimated$D <- mmm

## @knitr simulate
set.seed(1)
sim_gp <- lapply(1:100, function(i) ForwardSimulate(f, p, x_grid, h_grid, K, opt_gp$D, z_g, profit=profit))
set.seed(1)
sim_true <- lapply(1:100, function(i) ForwardSimulate(f, p, x_grid, h_grid, K, opt_true$D, z_g, profit=profit))
set.seed(1)
sim_est <- lapply(1:100, function(i) ForwardSimulate(f, p, x_grid, h_grid, K, opt_estimated$D, z_g, profit=profit))

## @knitr tidy
dat <- list(GP = sim_gp, Parametric = sim_est, True = sim_true)
dat <- melt(dat, id=names(dat[[1]][[1]]))
dt <- data.table(dat)
setnames(dt, c("L1", "L2"), c("method", "reps"))

## @knitr sim-fish
ggplot(dt) +
  geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=.1) +
  scale_colour_manual(values=cbPalette, guide = guide_legend(override.aes = list(alpha = 1)))

## @knitr sim-harvest
ggplot(dt) +
  geom_line(aes(time, harvest, group=interaction(reps,method), color=method), alpha=.1) +
  scale_colour_manual(values=cbPalette, guide = guide_legend(override.aes = list(alpha = 1)))

## @knitr costs
profits <- dt[, sum(profit), by = c("reps", "method")]
means <- profits[, mean(V1), by = method]
sds <- profits[, sd(V1), by = method]
yield <- cbind(means, sd = sds$V1)
yield


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

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

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

```{r May}
```

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

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

With parameters `r p`.

```{r x0-obs}
xT <- x_grid[2]
x_0_observed <- x_grid[60]
```

```{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-fit, dependson="lag-data"}
```

We fit a Gaussian process with

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

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

## 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 simulate}
```

```{r tidy}
```

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

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

```{r costs}
```

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


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_knit$set(fig.path="figure/small-noise/")
opts_chunk$set(fig.width=6, fig.height=6, echo=FALSE, dev="pdf", warning=FALSE, message=FALSE, 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
}

load("small-noise.rda")

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

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

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

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

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

save(list=ls(), file='small-noise.rda')


text:
```{r include=FALSE}
require(MASS)
require(ggplot2)
require(kernlab)
require(nonparametricbayes)
opts_knit$set(upload.fun = socialR::notebook.url)
opts_chunk$set(dev.args=list(bg="transparent"), comment=NA, tidy=FALSE, message=FALSE)
theme_update(plot.background = element_rect(fill = "transparent",colour = NA))
```

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

```{r results="hide", eval=FALSE, include=FALSE}
#knit("beverton_holt_data.Rmd")
#true <- data.frame(x = X, y = sapply(X, f, 0, p))
```

### Cholesky method

```{r}
  SE <- function(Xi,Xj, l=1) 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)
  L <- chol(K + sigma_n^2 * I)
  alpha <- forwardsolve(t(L), backsolve(L, obs$y)) # see http://stat.ethz.ch/R-manual/R-patched/library/base/html/chol2inv.html
  k_star <- cov(obs$x, X)
  Y <- t(k_star) %*% alpha
  v <- backsolve(L, k_star)
  Var <- cov(X,X) - t(v) %*% v
  loglik <- -.5 * t(obs$y) %*% alpha - sum(log(diag(L))) - n * log(2 * pi) / 2
```

### Direct method

``` {r}
  d <- gp_fit(obs, X, c(sigma_n=sigma_n, l=l, tau=1), "direct")
  Ef <- d$Ef
  Cf <- d$Cf
```

### Direct sequential method, avoids matrix inverse instability

``` {r}
s <- gp_fit(obs, X, c(sigma_n=sigma_n, l=l, tau=1), "sequential")
ef <- s$Ef
cf <- s$Cf

```

``` {r}
ggplot(data.frame(x=X, Ef=Ef, ef=ef))+ geom_point(aes(x,Ef), col='red') + geom_line(aes(x,ef))
````

### kernlab method

``` {r}
k <- gp_fit(obs, X, c(sigma_n=sigma_n, l=l, tau=1), "kernlab")
y_p <- k$Ef
```

Compare these results:

``` {r}
require(reshape2)
df <- data.frame(x = X, direct = Ef, Cholesky = Y, kernlab = y_p, sequential = ef)
df <- melt(df, id = "x")
ggplot(df)+ geom_jitter(aes(x, value, color = variable)) + geom_point(data = obs, aes(x,y))
````


text:
scenarios <- melt(parameters)

ids <- c(rep(1:3, 3), rep(4:6, 3))

scenarios <- cbind("scenario"= ids, scenarios)
names(scenarios) <- c("scenario", "value", "parameter", "model")
scenarios <- scenarios[c("scenario", "model", "parameter", "value")]

Allen <- scenarios$model == "Allen"
Myers <- scenarios$model == "Myers"

allen_ids <- scenarios[Allen, "parameter"]
myers_ids <- scenarios[Myers, "parameter"]

allen_names <- names(parameters$Allen[[1]])
myers_names <- names(parameters$Myers[[1]])

scenarios[Allen, "parameter"] <- allen_names[allen_ids]
scenarios[Myers, "parameter"] <- myers_names[myers_ids]

scenarios


text:
```{r, include=FALSE}
## @knitr plotting-options
source("~/.knitr_defaults.R")
opts_knit$set(upload.fun = socialR::flickr.url)
library(knitcitations)
library(nonparametricbayes)
opts_chunk$set(cache=FALSE)
```

```{r}
f <- RickerAllee
p <- c(.5, .8, .5)
K <- .8  # approx, a li'l' less
allee <- .5 # approx, a li'l' less
```

```{r}
sigma_g <- 0.005
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+.2# observations start from
x0 <- K # simulation under policy starts from
Tobs <- 40
MaxT = 1000 # timeout for value iteration convergence
```

```{r}
set.seed(1234)
x <- numeric(Tobs)
x[1] <- Xo
nz <- 1
for(t in 1:(Tobs-1))
  x[t+1] = z_g() * f(x[t], h=0, p=p)
obs <- data.frame(x = c(rep(0,nz),
                        pmax(rep(0,Tobs-1), x[1:(Tobs-1)])),
                  y = c(rep(0,nz),
                        x[2:Tobs]))
raw_plot <- ggplot(data.frame(time = 1:Tobs, x=x), aes(time,x)) + geom_line()
raw_plot
```

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

## @knitr mle-output
true_means <- sapply(x_grid, f, 0, p)
est_means <- sapply(x_grid, est$f, 0, est$p)

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

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

```{r}
## @knitr gp_traces_densities
gp_assessment_plots <- summary_gp_mcmc(gp, burnin=1e4, thin=300)
gp_assessment_plots
```

```{r}
## @knitr gp-output
tgp_dat <- data.frame(x = x_grid,
               y = gp_dat$E_Ef,
               ymin = pmax(gp_dat$E_Ef - 2 * sqrt(gp_dat$E_Vf),0),
               ymax = gp_dat$E_Ef + 2 * sqrt(gp_dat$E_Vf) )

## @knitr assemble-models
models <- data.frame(x=x_grid,
                     GP=tgp_dat$y,
                     True=true_means)
models <- melt(models, id="x")
names(models) <- c("x", "method", "value")

## @knitr Figure1
plot_gp <- ggplot(tgp_dat) + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") +
  geom_line(data=models, aes(x, value, col=method), lwd=1, alpha=0.8) +
  geom_point(data=obs, aes(x,y), alpha=0.8) +
  xlab(expression(X[t])) + ylab(expression(X[t+1]))
print(plot_gp)
```

```{r}
dat_fig1c <- melt(lapply(sample(100),
  function(i){
    data.frame(stock = x_grid,
                GP = pmax( # non-negative curves only
                          mvrnorm(1,
                                  gp_dat$Ef_posterior[,i],
                                  gp_dat$E_Cf # could draw from Cf_posterior after assuring is pos-def
                                  )
                          ,0)
               )

  }), id="stock")

ggplot(dat_fig1c, aes(stock, value, group=L1)) + geom_line(alpha=.3)

# a <- dcast(dat_fig1c, stock ~ L1, value.var='value')
# dat_fig1c <- data.table(dat_fig1c)
# to_grid <- function(x, grid) sapply(x, function(x) which.min(abs(grid - x)))
# dat_fig1c$stock <- to_grid(dat_fig1c$stock, x_grid)
# a <- dat_fig1c[, to_grid(value, x_grid), by=c('stock')]
```

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

opt_gp <- find_dp_optim(matrices_gp, x_grid, h_grid, 5, xT=0.5, profit, delta, reward=1)
opt_gp$D <- opt_gp$D[,1]

## @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=.01)
```

```{r}
OPT = data.frame(GP = opt_gp$D, True = opt_true$D)

## @knitr Figure2
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")
```

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

## @knitr Figure3
ggplot(sims_data) +
  geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=.1)
```

```{r}
## @knitr profits
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]])
actual_over_optimal <-subset(tmp, variable != "True")

## @knitr Figure4
ggplot(actual_over_optimal, aes(value)) + geom_histogram(aes(fill=variable)) +
  facet_wrap(~variable, scales = "free_y")  + guides(legend.position = "none") +
  xlab("Total profit by replicate")
```

# Forecasts

## One-step ahead predictor plots and forecasting plots

Now we show how each of the fitted models performs on the training data (e.g. plot of the step-ahead predictors). For the GP, we need to predict on the training data first:

```{r}
gp_f_at_obs <- gp_predict(gp, x, burnin=1e4, thin=300)
```

For the parametric models, prediction is a matter of feeding in each of the observed data points to the fitted parameters, like so

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

Which we apply over each of the fitted models, including the GP, organizing the "expected" transition points (given the previous point) into a data frame.

```{r}
df <- melt(data.frame(time = 1:length(x), stock = x,
                GP = gp_f_at_obs$E_Ef,
                True = step_ahead(x,f,p)
                 ), id=c("time", "stock"))

ggplot(df) + geom_point(aes(time, stock)) +
  geom_line(aes(time, value, col=variable))
```

## Posterior predictive curves

This shows only the mean predictions.  For the Bayesian cases, we can instead loop over the posteriors of the parameters (or samples from the GP posterior) to get the distribution of such curves in each case.

We will need a vector version (`pmin` in place of `min`) of this function that can operate on the posteriors, others are vectorized already.

```{r}
ricker_f <- function(x,h,p){
  sapply(x, function(x){
    x <- pmax(0, x-h)
    pmax(0, x * exp(p[2] * (1 - x / p[1] )) )
  })
}
```

Then we proceed as before, now looping over 100 random samples from the posterior for each Bayesian estimate.  We write this as a function for easy reuse.

```{r}
require(MASS)

step_ahead_posteriors <- function(x){
gp_f_at_obs <- gp_predict(gp, x, burnin=1e4, thin=300)
df_post <- melt(lapply(sample(100),
  function(i){
    data.frame(time = 1:length(x), stock = x,
                GP = pmax(mvrnorm(1,
                                  gp_f_at_obs$Ef_posterior[,i],
                                  pmax(gp_f_at_obs$E_Cf, 1e-4)
                                  #pmax(gp_f_at_obs$Cf_posterior[[i]],1e-5)
                                  ),0),
                True = step_ahead(x,f,p))
  }), id=c("time", "stock"))

}

df_post <- step_ahead_posteriors(x)

ggplot(df_post) + geom_point(aes(time, stock)) +
  geom_line(aes(time, value, col=variable, group=interaction(L1,variable)), alpha=.1)
```

## Performance on data outside of observations

Of course it is hardly suprising that all models do reasonably well on the data on which they were trained.  A crux of the problem is the model performance on data outside the observed range.  (Though we might also wish to repeat the above plot on data in the observed range but a different sequence from the observed data).

First we generate some data from the underlying model coming from below the tipping point:

```{r}
Tobs <- 30
y <- numeric(Tobs)
y[1] = .45
for(t in 1:(Tobs-1))
    y[t+1] = z_g() * f(y[t], h=0, p=p)
```

Proceed as before on this data:

```{r}
crash_data <- step_ahead_posteriors(y)

ggplot(crash_data) + geom_point(aes(time, stock)) +
  geom_line(aes(time, value, col=variable, group=interaction(L1,variable)), alpha=.1) +
  facet_wrap(~variable)

```

## Forecast Distributions

Another way to visualize this is to look directly at the distribution predicted under each model, one step and several steps into the future (say, at a fixed harvest level).

```{r}
library(expm)
get_forecasts <- function(i, Tobs, h_i){

df <- data.frame(
  x = x_grid,
  GP = matrices_gp[[h_i]][i,],
  True = matrices_true[[h_i]][i,])

df2 <- data.frame(
  x = x_grid,
  GP = (matrices_gp[[h_i]] %^% Tobs)[i,],
  True = (matrices_true[[h_i]] %^% Tobs)[i,])

T2 <- 2 * Tobs

df4 <- data.frame(
  x = x_grid,
  GP = (matrices_gp[[h_i]] %^% T2)[i,],
  True = (matrices_true[[h_i]] %^% T2)[i,])

forecasts <- melt(list(T_start = df, T_mid = df2, T_end = df4), id="x")

}

i = 15
```

This takes `i`, an index to `x_grid` value (e.g. for `i`=`r i` corresponds to a starting postion `x` = `r x_grid[i]`)

```{r}
forecasts <- get_forecasts(i = 15, Tobs = 15, h_i = 1)
ggplot(forecasts) +
  geom_line(aes(x, value, group=interaction(variable, L1), col=variable, lty=L1)) +
  facet_wrap(~ variable, scale="free_y", ncol=2)
```

We can compare to a better starting stock,

```{r}
i<-45
```

where `i`=`r i` corresponds to a starting postion `x` = `r x_grid[i]`

```{r}
forecasts <- get_forecasts(i = i, Tobs = 10, h_i = 1)

ggplot(forecasts) +
  geom_line(aes(x, value, group=interaction(variable, L1), col=variable, lty=L1)) +
  facet_wrap(~ variable, scale="free_y", ncol=2)
```

and a middle value showing the population decline:

```{r}
i<-35
```

where `i`=`r i` corresponds to a starting postion `x` = `r x_grid[i]`

```{r}
forecasts <- get_forecasts(i = i, Tobs = 10, h_i = 1)

ggplot(forecasts) +
  geom_line(aes(x, value, group=interaction(variable, L1), col=variable, lty=L1)) +
  facet_wrap(~ variable, scale="free_y", ncol=2)
```


text:
```r
f <- RickerAllee
p <- c(.5, .8, .5)
K <- .8  # approx, a li'l' less
allee <- .5 # approx, a li'l' less
```

```r
sigma_g <- 0.005
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+.2# observations start from
x0 <- K # simulation under policy starts from
Tobs <- 40
MaxT = 1000 # timeout for value iteration convergence
```

```r
set.seed(1234)
x <- numeric(Tobs)
x[1] <- Xo
nz <- 1
for(t in 1:(Tobs-1))
  x[t+1] = z_g() * f(x[t], h=0, p=p)
obs <- data.frame(x = c(rep(0,nz),
                        pmax(rep(0,Tobs-1), x[1:(Tobs-1)])),
                  y = c(rep(0,nz),
                        x[2:Tobs]))
raw_plot <- ggplot(data.frame(time = 1:Tobs, x=x), aes(time,x)) + geom_line()
raw_plot
```

![plot of chunk unnamed-chunk-4](http://farm6.staticflickr.com/5513/9488645423_86e44c3679_o.png)

```r
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)

## @knitr mle-output
true_means <- sapply(x_grid, f, 0, p)
est_means <- sapply(x_grid, est$f, 0, est$p)

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

```r
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)
```

```r
## @knitr gp_traces_densities
gp_assessment_plots <- summary_gp_mcmc(gp, burnin=1e4, thin=300)
gp_assessment_plots
```

```
$traces_plot
```

![plot of chunk unnamed-chunk-7](http://farm4.staticflickr.com/3802/9491452730_34ef90a140_o.png)

```

$posteriors_plot
```

![plot of chunk unnamed-chunk-7](http://farm3.staticflickr.com/2811/9488657493_ffd898069f_o.png)

```r
## @knitr gp-output
tgp_dat <- data.frame(x = x_grid,
               y = gp_dat$E_Ef,
               ymin = pmax(gp_dat$E_Ef - 2 * sqrt(gp_dat$E_Vf),0),
               ymax = gp_dat$E_Ef + 2 * sqrt(gp_dat$E_Vf) )

## @knitr assemble-models
models <- data.frame(x=x_grid,
                     GP=tgp_dat$y,
                     True=true_means)
models <- melt(models, id="x")
names(models) <- c("x", "method", "value")

## @knitr Figure1
plot_gp <- ggplot(tgp_dat) + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") +
  geom_line(data=models, aes(x, value, col=method), lwd=1, alpha=0.8) +
  geom_point(data=obs, aes(x,y), alpha=0.8) +
  xlab(expression(X[t])) + ylab(expression(X[t+1]))
print(plot_gp)
```

![plot of chunk unnamed-chunk-8](http://farm4.staticflickr.com/3681/9488658009_efdb474fe9_o.png)

```r
dat_fig1c <- melt(lapply(sample(100),
  function(i){
    data.frame(stock = x_grid,
                GP = pmax( # non-negative curves only
                          mvrnorm(1,
                                  gp_dat$Ef_posterior[,i],
                                  gp_dat$E_Cf # could draw from Cf_posterior after assuring is pos-def
                                  )
                          ,0)
               )

  }), id="stock")

ggplot(dat_fig1c, aes(stock, value, group=L1)) + geom_line(alpha=.3)
```

![plot of chunk unnamed-chunk-9](http://farm3.staticflickr.com/2848/9491454504_4b8d995bd3_o.png)

```r

# a <- dcast(dat_fig1c, stock ~ L1, value.var='value')
# dat_fig1c <- data.table(dat_fig1c)
# to_grid <- function(x, grid) sapply(x, function(x) which.min(abs(grid - x)))
# dat_fig1c$stock <- to_grid(dat_fig1c$stock, x_grid)
# a <- dat_fig1c[, to_grid(value, x_grid), by=c('stock')]
```

```r
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)

opt_gp <- find_dp_optim(matrices_gp, x_grid, h_grid, 5, xT=0.5, profit, delta, reward=1)
opt_gp$D <- opt_gp$D[,1]

## @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=.01)
```

```r
OPT = data.frame(GP = opt_gp$D, True = opt_true$D)

## @knitr Figure2
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")
```

![plot of chunk unnamed-chunk-11](http://farm4.staticflickr.com/3683/9488670189_2e7168e731_o.png)

```r
## @knitr sims
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))

## @knitr Figure3
ggplot(sims_data) +
  geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=.1)
```

![plot of chunk unnamed-chunk-12](http://farm4.staticflickr.com/3728/9488671349_5612a1de51_o.png)

```r
## @knitr profits
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]])
actual_over_optimal <-subset(tmp, variable != "True")

## @knitr Figure4
ggplot(actual_over_optimal, aes(value)) + geom_histogram(aes(fill=variable)) +
  facet_wrap(~variable, scales = "free_y")  + guides(legend.position = "none") +
  xlab("Total profit by replicate")
```

![plot of chunk unnamed-chunk-13](http://farm4.staticflickr.com/3761/9488672199_2355ecaa47_o.png)

# Forecasts

## One-step ahead predictor plots and forecasting plots

Now we show how each of the fitted models performs on the training data (e.g. plot of the step-ahead predictors). For the GP, we need to predict on the training data first:

```r
gp_f_at_obs <- gp_predict(gp, x, burnin=1e4, thin=300)
```

For the parametric models, prediction is a matter of feeding in each of the observed data points to the fitted parameters, like so

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

Which we apply over each of the fitted models, including the GP, organizing the "expected" transition points (given the previous point) into a data frame.

```r
df <- melt(data.frame(time = 1:length(x), stock = x,
                GP = gp_f_at_obs$E_Ef,
                True = step_ahead(x,f,p)
                 ), id=c("time", "stock"))

ggplot(df) + geom_point(aes(time, stock)) +
  geom_line(aes(time, value, col=variable))
```

![plot of chunk unnamed-chunk-16](http://farm8.staticflickr.com/7407/9488672745_58e2364a6a_o.png)

## Posterior predictive curves

This shows only the mean predictions.  For the Bayesian cases, we can instead loop over the posteriors of the parameters (or samples from the GP posterior) to get the distribution of such curves in each case.

We will need a vector version (`pmin` in place of `min`) of this function that can operate on the posteriors, others are vectorized already.

```r
ricker_f <- function(x,h,p){
  sapply(x, function(x){
    x <- pmax(0, x-h)
    pmax(0, x * exp(p[2] * (1 - x / p[1] )) )
  })
}
```

Then we proceed as before, now looping over 100 random samples from the posterior for each Bayesian estimate.  We write this as a function for easy reuse.

```r
require(MASS)

step_ahead_posteriors <- function(x){
gp_f_at_obs <- gp_predict(gp, x, burnin=1e4, thin=300)
df_post <- melt(lapply(sample(100),
  function(i){
    data.frame(time = 1:length(x), stock = x,
                GP = pmax(mvrnorm(1,
                                  gp_f_at_obs$Ef_posterior[,i],
                                  pmax(gp_f_at_obs$E_Cf, 1e-4)
                                  #pmax(gp_f_at_obs$Cf_posterior[[i]],1e-5)
                                  ),0),
                True = step_ahead(x,f,p))
  }), id=c("time", "stock"))

}

df_post <- step_ahead_posteriors(x)

ggplot(df_post) + geom_point(aes(time, stock)) +
  geom_line(aes(time, value, col=variable, group=interaction(L1,variable)), alpha=.1)
```

![plot of chunk unnamed-chunk-18](http://farm3.staticflickr.com/2828/9491469452_02e30cee68_o.png)

## Performance on data outside of observations

Of course it is hardly suprising that all models do reasonably well on the data on which they were trained.  A crux of the problem is the model performance on data outside the observed range.  (Though we might also wish to repeat the above plot on data in the observed range but a different sequence from the observed data).

First we generate some data from the underlying model coming from below the tipping point:

```r
Tobs <- 30
y <- numeric(Tobs)
y[1] = .45
for(t in 1:(Tobs-1))
    y[t+1] = z_g() * f(y[t], h=0, p=p)
```

Proceed as before on this data:

```r
crash_data <- step_ahead_posteriors(y)
```

```
Error: 'Sigma' is not positive definite
```

```r

ggplot(crash_data) + geom_point(aes(time, stock)) +
  geom_line(aes(time, value, col=variable, group=interaction(L1,variable)), alpha=.1) +
  facet_wrap(~variable)
```

```
Error: object 'crash_data' not found
```

## Forecast Distributions

Another way to visualize this is to look directly at the distribution predicted under each model, one step and several steps into the future (say, at a fixed harvest level).

```r
library(expm)
get_forecasts <- function(i, Tobs, h_i){

df <- data.frame(
  x = x_grid,
  GP = matrices_gp[[h_i]][i,],
  True = matrices_true[[h_i]][i,])

df2 <- data.frame(
  x = x_grid,
  GP = (matrices_gp[[h_i]] %^% Tobs)[i,],
  True = (matrices_true[[h_i]] %^% Tobs)[i,])

T2 <- 2 * Tobs

df4 <- data.frame(
  x = x_grid,
  GP = (matrices_gp[[h_i]] %^% T2)[i,],
  True = (matrices_true[[h_i]] %^% T2)[i,])

forecasts <- melt(list(T_start = df, T_mid = df2, T_end = df4), id="x")

}

i = 15
```

This takes `i`, an index to `x_grid` value (e.g. for `i`=15 corresponds to a starting postion `x` = 0.3429)

```r
forecasts <- get_forecasts(i = 15, Tobs = 15, h_i = 1)
ggplot(forecasts) +
  geom_line(aes(x, value, group=interaction(variable, L1), col=variable, lty=L1)) +
  facet_wrap(~ variable, scale="free_y", ncol=2)
```

![plot of chunk unnamed-chunk-22](http://farm3.staticflickr.com/2834/9488674537_253891e90f_o.png)

We can compare to a better starting stock,

```r
i<-45
```

where `i`=45 corresponds to a starting postion `x` = 1.0776

```r
forecasts <- get_forecasts(i = i, Tobs = 10, h_i = 1)

ggplot(forecasts) +
  geom_line(aes(x, value, group=interaction(variable, L1), col=variable, lty=L1)) +
  facet_wrap(~ variable, scale="free_y", ncol=2)
```

![plot of chunk unnamed-chunk-24](http://farm6.staticflickr.com/5519/9488674991_bd7fd57cfc_o.png)

and a middle value showing the population decline:

```r
i<-35
```

where `i`=35 corresponds to a starting postion `x` = 0.8327

```r
forecasts <- get_forecasts(i = i, Tobs = 10, h_i = 1)

ggplot(forecasts) +
  geom_line(aes(x, value, group=interaction(variable, L1), col=variable, lty=L1)) +
  facet_wrap(~ variable, scale="free_y", ncol=2)
```

![plot of chunk unnamed-chunk-26](http://farm8.staticflickr.com/7290/9488675377_3b9fb86078_o.png)


text:
Fixed priors on hyperparameters, fixed model type.

```r
s2.p <- c(50,50)  # inverse gamma mean is
tau2.p <- c(20,1)
d.p = c(10, 1/0.01, 10, 1/0.01)
nug.p = c(10, 1/0.01, 10, 1/0.01) # 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
z_g = function() rlnorm(1, 0, sigma_g)
z_m = function() 1+(2*runif(1, 0,  1)-1) * sigma_m
```

```r
f <- RickerAllee
p <- c(1.1, 10, 5)
K <- 10
allee <- 5
# valid only for p[2] = 2
#K <- p[1] * p[3] / 2 + sqrt( (p[1] * p[3]) ^ 2 - 4 * p[3] ) / 2 #
#allee <- p[1] * p[3] / 2 - sqrt( (p[1] * p[3]) ^ 2 - 4 * p[3] ) / 2 # allee threshold
#e_star <- (p[1] * sqrt(p[3]) - 2) / 2 ## Bifurcation point
```

```r
sigma_g <- 0.05
sigma_m <- 0.02
x_grid <- seq(0, 1.5 * K, length=101)
h_grid <- x_grid
```

With parameters 1.1, 10, 5.

```r

sim_obs <- function(Xo, z_g, f, p, Tobs = 40, seed = 1){
  x <- numeric(Tobs)
  x[1] <- Xo
  set.seed(seed)
  for(t in 1:(Tobs-1))
    x[t+1] = z_g() * f(x[t], h=0, p=p)
  data_plot <- qplot(1:Tobs, x) + ggtitle(paste("seed", seed))
  print(data_plot)
  obs <- data.frame(x=c(0,x[1:(Tobs-1)]),y=c(0,x[2:Tobs]))
  obs
}

par_est <- function(obs){
  estf <- function(p){
    mu <- log(obs$x) + p["r"]*(1-obs$x/p["K"])
    llik <- -sum(dlnorm(obs$y, mu, p["s"]), log=TRUE)
    if(!is.numeric(llik) | is.nan(llik) | !(llik < Inf)){
      warning("possible error in llik")
      llik <- 1e30
    }
    llik
  }
  o <- optim(par = c(r=1,K=mean(obs$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']
  list(f_alt = f_alt, p_alt = p_alt, sigma_g_alt = sigma_g_alt)
}

optimal_policy <- function(gp, f, f_alt, p, p_alt, x_grid, h_grid, sigma_g, sigma_g_alt, delta, xT, profit, reward, OptTime){
  matrices_gp <- gp_transition_matrix(gp$ZZ.km, gp$ZZ.ks2, x_grid, h_grid)
  opt_gp <- find_dp_optim(matrices_gp, x_grid, h_grid, OptTime, xT, profit, delta, reward=reward)
  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)
  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)
  list(gp_D = opt_gp$D[,1], true_D = opt_true$D[,1], est_D = opt_estimated$D[,1])
}

simulate_opt <- function(OPT, f, p, x_grid, h_grid, x0, z_g, profit){
  gp_D <- sapply(1:OptTime, function(i) OPT$gp_D)
  true_D <- sapply(1:OptTime, function(i) OPT$true_D)
  est_D <- sapply(1:OptTime, function(i) OPT$est_D)
  set.seed(1)
  sim_gp <- lapply(1:100, function(i) ForwardSimulate(f, p, x_grid, h_grid, x0, gp_D, z_g, profit=profit))
  set.seed(1)
  sim_true <- lapply(1:100, function(i) ForwardSimulate(f, p, x_grid, h_grid, x0, true_D, z_g, profit=profit))
  set.seed(1)
  sim_est <- lapply(1:100, function(i) ForwardSimulate(f, p, x_grid, h_grid, x0, est_D, z_g, profit=profit))
  dat <- list(GP = sim_gp, Parametric = sim_est, True = sim_true)
  dat <- melt(dat, id=names(dat[[1]][[1]]))
  dt <- data.table(dat)
  setnames(dt, c("L1", "L2"), c("method", "reps"))
  dt
}

sim_plots <- function(dt, seed=1){
  fish_plot <- ggplot(dt) +
    geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=.1) +
    scale_colour_manual(values=cbPalette, guide = guide_legend(override.aes = list(alpha = 1)))+
    ggtitle(paste("seed", seed))
  print(fish_plot)
  harvest_plot <- ggplot(dt) +
    geom_line(aes(time, harvest, group=interaction(reps,method), color=method), alpha=.1) +
    scale_colour_manual(values=cbPalette, guide = guide_legend(override.aes = list(alpha = 1))) +
    ggtitle(paste("seed", seed))
  print(harvest_plot)
}

profits_stats <- function(dt){
  profits <- dt[, sum(profit), by = c("reps", "method")]
  means <- profits[, mean(V1), by = method]
  sds <- profits[, sd(V1), by = method]
  yield <- cbind(means, sd = sds$V1)
  yield
}

gp_plot <- function(gp, f, p, f_alt, p_alt, x_grid, obs, seed){
  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))
  true <- sapply(x_grid, f, 0, p)
  est <- sapply(x_grid, f_alt, 0, p_alt)
  models <- data.frame(x=x_grid, GP=tgp_dat$y, Parametric=est, True=true)
  models <- melt(models, id="x")
  names(models) <- c("x", "method", "value")
  plot_gp <- ggplot(tgp_dat) + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") +
    geom_line(data=models, aes(x, value, col=method), lwd=2, alpha=0.8) +
    geom_point(data=obs, aes(x,y), alpha=0.8) +
    xlab(expression(X[t])) + ylab(expression(X[t+1])) +
    scale_colour_manual(values=cbPalette) +
    ggtitle(paste("seed", seed))
  print(plot_gp)
}

posteriors_plot <- function(gp, priors){
  hyperparameters <- c("index", "s2", "tau2", "beta0", "nug", "d", "ldetK")
  posteriors <- melt(gp$trace$XX[[1]][,hyperparameters], id="index")
  prior_curves <- ddply(posteriors, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = priors[[dd$variable[1]]](grid))
  })
  plot_posteriors <- ggplot(posteriors) +
    #geom_density(aes(value), lwd=2) +
    geom_histogram(aes(x=value, y=..density..), alpha=0.7) +
    geom_line(data=prior_curves, aes(x=value, y=density), col="red", lwd=2) +
    facet_wrap(~ variable, scale="free")
  print(plot_posteriors)
}
```

```r
require(snowfall)
sfInit(par=TRUE, cpu=8)
sfExportAll()
sfLibrary(pdgControl)
sfLibrary(nonparametricbayes)
sfLibrary(reshape2)
sfLibrary(ggplot2)
sfLibrary(data.table)
sfLibrary(tgp)
sfLibrary(kernlab)
sfLibrary(MCMCpack)
sfLibrary(plyr)
```

```r
seed <- round(runif(32) * 1e6)
seed[1] <- 1
seed
```

```
 [1]      1 969343 449127 263444 674131 949894 913719 444907 989690 229414
[11] 865064 804050 965462 777512 497134  65605 141485 193569 259734 451474
[21] 695013 530379 124241 497581 611219 718556 620409  26555  71806 717469
[31] 393034 540837
```

```r
yields <- lapply(seed,
function(seed_i){
  Xo <- allee + x_grid[10] # observations start from
  x0 <- Xo # simulation under policy starts from
  obs <- sim_obs(Xo, z_g, f, p, Tobs=40, seed = seed_i)
  est <- par_est(obs)
  gp <- bgp(X=obs$x, XX=x_grid, Z=obs$y, verb=0,
          meanfn="constant", bprior="b0", BTE=c(2000,16000,2),
          m0r1=FALSE, corr="exp", trace=FALSE,
          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_alt, est$p_alt, x_grid, obs, seed_i)
#  posteriors_plot(gp, priors) # needs trace=TRUE!
  OPT <- optimal_policy(gp, f, est$f_alt, p, est$p_alt, x_grid, h_grid, sigma_g, est$sigma_g_alt, delta, xT, profit, reward, OptTime)
  dt <- simulate_opt(OPT, f, p, x_grid, h_grid, x0, z_g, profit)
  sim_plots(dt, seed=seed_i)
  profits_stats(dt)
})
```

![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-06-60abafd424-unnamed-chunk-21.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-06-60abafd424-unnamed-chunk-22.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-07-60abafd424-unnamed-chunk-23.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-08-60abafd424-unnamed-chunk-24.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-09-60abafd424-unnamed-chunk-25.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-09-60abafd424-unnamed-chunk-26.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-10-60abafd424-unnamed-chunk-27.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-11-60abafd424-unnamed-chunk-28.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-12-60abafd424-unnamed-chunk-29.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-12-60abafd424-unnamed-chunk-210.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-13-60abafd424-unnamed-chunk-211.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-14-60abafd424-unnamed-chunk-212.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-15-60abafd424-unnamed-chunk-213.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-15-60abafd424-unnamed-chunk-214.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-16-60abafd424-unnamed-chunk-215.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-17-60abafd424-unnamed-chunk-216.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-18-60abafd424-unnamed-chunk-217.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-18-60abafd424-unnamed-chunk-218.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-19-60abafd424-unnamed-chunk-219.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-20-60abafd424-unnamed-chunk-220.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-21-60abafd424-unnamed-chunk-221.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-21-60abafd424-unnamed-chunk-222.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-22-60abafd424-unnamed-chunk-223.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-23-60abafd424-unnamed-chunk-224.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-23-60abafd424-unnamed-chunk-225.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-24-60abafd424-unnamed-chunk-226.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-25-60abafd424-unnamed-chunk-227.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-26-60abafd424-unnamed-chunk-228.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-26-60abafd424-unnamed-chunk-229.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-27-60abafd424-unnamed-chunk-230.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-28-60abafd424-unnamed-chunk-231.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-29-60abafd424-unnamed-chunk-232.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-29-60abafd424-unnamed-chunk-233.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-30-60abafd424-unnamed-chunk-234.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-31-60abafd424-unnamed-chunk-235.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-32-60abafd424-unnamed-chunk-236.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-32-60abafd424-unnamed-chunk-237.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-33-60abafd424-unnamed-chunk-238.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-34-60abafd424-unnamed-chunk-239.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-35-60abafd424-unnamed-chunk-240.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-35-60abafd424-unnamed-chunk-241.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-36-60abafd424-unnamed-chunk-242.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-37-60abafd424-unnamed-chunk-243.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-37-60abafd424-unnamed-chunk-244.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-38-60abafd424-unnamed-chunk-245.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-39-60abafd424-unnamed-chunk-246.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-40-60abafd424-unnamed-chunk-247.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-40-60abafd424-unnamed-chunk-248.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-41-60abafd424-unnamed-chunk-249.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-42-60abafd424-unnamed-chunk-250.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-43-60abafd424-unnamed-chunk-251.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-43-60abafd424-unnamed-chunk-252.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-44-60abafd424-unnamed-chunk-253.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-45-60abafd424-unnamed-chunk-254.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-45-60abafd424-unnamed-chunk-255.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-46-60abafd424-unnamed-chunk-256.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-47-60abafd424-unnamed-chunk-257.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-53-60abafd424-unnamed-chunk-258.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-54-60abafd424-unnamed-chunk-259.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-54-60abafd424-unnamed-chunk-260.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-55-60abafd424-unnamed-chunk-261.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-56-60abafd424-unnamed-chunk-262.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-56-60abafd424-unnamed-chunk-263.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-57-60abafd424-unnamed-chunk-264.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-58-60abafd424-unnamed-chunk-265.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-59-60abafd424-unnamed-chunk-266.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-12-59-60abafd424-unnamed-chunk-267.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-00-60abafd424-unnamed-chunk-268.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-01-60abafd424-unnamed-chunk-269.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-02-60abafd424-unnamed-chunk-270.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-02-60abafd424-unnamed-chunk-271.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-03-60abafd424-unnamed-chunk-272.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-04-60abafd424-unnamed-chunk-273.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-05-60abafd424-unnamed-chunk-274.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-05-60abafd424-unnamed-chunk-275.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-06-60abafd424-unnamed-chunk-276.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-07-60abafd424-unnamed-chunk-277.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-08-60abafd424-unnamed-chunk-278.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-08-60abafd424-unnamed-chunk-279.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-09-60abafd424-unnamed-chunk-280.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-10-60abafd424-unnamed-chunk-281.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-11-60abafd424-unnamed-chunk-282.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-11-60abafd424-unnamed-chunk-283.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-12-60abafd424-unnamed-chunk-284.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-13-60abafd424-unnamed-chunk-285.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-14-60abafd424-unnamed-chunk-286.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-14-60abafd424-unnamed-chunk-287.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-15-60abafd424-unnamed-chunk-288.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-16-60abafd424-unnamed-chunk-289.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-17-60abafd424-unnamed-chunk-290.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-17-60abafd424-unnamed-chunk-291.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-18-60abafd424-unnamed-chunk-292.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-19-60abafd424-unnamed-chunk-293.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-20-60abafd424-unnamed-chunk-294.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-20-60abafd424-unnamed-chunk-295.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-21-60abafd424-unnamed-chunk-296.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-22-60abafd424-unnamed-chunk-297.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-23-60abafd424-unnamed-chunk-298.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-23-60abafd424-unnamed-chunk-299.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-24-60abafd424-unnamed-chunk-2100.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-25-60abafd424-unnamed-chunk-2101.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-26-60abafd424-unnamed-chunk-2102.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-26-60abafd424-unnamed-chunk-2103.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-27-60abafd424-unnamed-chunk-2104.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-28-60abafd424-unnamed-chunk-2105.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-29-60abafd424-unnamed-chunk-2106.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-29-60abafd424-unnamed-chunk-2107.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-30-60abafd424-unnamed-chunk-2108.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-31-60abafd424-unnamed-chunk-2109.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-31-60abafd424-unnamed-chunk-2110.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-32-60abafd424-unnamed-chunk-2111.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-33-60abafd424-unnamed-chunk-2112.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-34-60abafd424-unnamed-chunk-2113.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-34-60abafd424-unnamed-chunk-2114.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-35-60abafd424-unnamed-chunk-2115.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-36-60abafd424-unnamed-chunk-2116.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-37-60abafd424-unnamed-chunk-2117.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-37-60abafd424-unnamed-chunk-2118.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-38-60abafd424-unnamed-chunk-2119.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-39-60abafd424-unnamed-chunk-2120.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-40-60abafd424-unnamed-chunk-2121.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-41-60abafd424-unnamed-chunk-2122.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-41-60abafd424-unnamed-chunk-2123.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-42-60abafd424-unnamed-chunk-2124.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-43-60abafd424-unnamed-chunk-2125.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-44-60abafd424-unnamed-chunk-2126.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-44-60abafd424-unnamed-chunk-2127.png) ![plot of chunk unnamed-chunk-2](http://carlboettiger.info/assets/figures/2012-12-20-23-13-45-60abafd424-unnamed-chunk-2128.png)

```r
yields <- melt(yields, id=c("method", "V1", "sd"))
ggplot(yields) + geom_histogram(aes(V1)) + facet_wrap(~method, scale='free')
```

```
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-3](http://carlboettiger.info/assets/figures/2012-12-20-23-13-46-60abafd424-unnamed-chunk-3.png)

```r
yields
```

```
       method     V1      sd L1
1          GP 8.1300 1.69491  1
2  Parametric 2.4240 0.27111  1
3        True 7.8840 1.76601  1
4          GP 8.0174 1.67504  2
5  Parametric 2.6360 0.34841  2
6        True 7.8840 1.76601  2
7          GP 6.4140 1.44192  3
8  Parametric 0.5685 0.36637  3
9        True 7.8840 1.76601  3
10         GP 7.6080 1.55418  4
11 Parametric 0.1485 0.15739  4
12       True 7.8840 1.76601  4
13         GP 8.0955 1.69390  5
14 Parametric 0.0915 0.11850  5
15       True 7.8840 1.76601  5
16         GP 6.3000 0.00000  6
17 Parametric 2.7797 0.22161  6
18       True 7.8840 1.76601  6
19         GP 8.0745 1.66938  7
20 Parametric 2.5102 0.41036  7
21       True 7.8840 1.76601  7
22         GP 6.3000 0.00000  8
23 Parametric 2.8819 1.13615  8
24       True 7.8840 1.76601  8
25         GP 8.0775 1.66520  9
26 Parametric 0.5685 0.36637  9
27       True 7.8840 1.76601  9
28         GP 8.1225 1.69992 10
29 Parametric 2.7786 0.08342 10
30       True 7.8840 1.76601 10
31         GP 8.0759 1.66078 11
32 Parametric 2.5931 0.36977 11
33       True 7.8840 1.76601 11
34         GP 8.0835 1.72560 12
35 Parametric 6.9945 1.95619 12
36       True 7.8840 1.76601 12
37         GP 7.3545 1.49506 13
38 Parametric 2.3198 0.41146 13
39       True 7.8840 1.76601 13
40         GP 8.0775 1.66520 14
41 Parametric 2.5945 0.23523 14
42       True 7.8840 1.76601 14
43         GP 6.3000 0.00000 15
44 Parametric 0.5685 0.36637 15
45       True 7.8840 1.76601 15
46         GP 6.9930 1.46548 16
47 Parametric 2.5945 0.23523 16
48       True 7.8840 1.76601 16
49         GP 8.1120 1.69715 17
50 Parametric 2.3114 0.23917 17
51       True 7.8840 1.76601 17
52         GP 7.3200 1.51222 18
53 Parametric 2.6059 0.27427 18
54       True 7.8840 1.76601 18
55         GP 8.1180 1.69388 19
56 Parametric 0.5685 0.36637 19
57       True 7.8840 1.76601 19
58         GP 8.0220 1.74766 20
59 Parametric 5.5725 1.70885 20
60       True 7.8840 1.76601 20
61         GP 7.6080 1.55418 21
62 Parametric 3.1926 1.10508 21
63       True 7.8840 1.76601 21
64         GP 6.4695 1.43474 22
65 Parametric 2.3919 0.39680 22
66       True 7.8840 1.76601 22
67         GP 7.9830 1.63145 23
68 Parametric 2.2469 0.16183 23
69       True 7.8840 1.76601 23
70         GP 8.1285 1.69568 24
71 Parametric 2.6167 0.22072 24
72       True 7.8840 1.76601 24
73         GP 6.9960 1.46792 25
74 Parametric 2.4249 0.12660 25
75       True 7.8840 1.76601 25
76         GP 6.3000 0.00000 26
77 Parametric 3.1950 0.95032 26
78       True 7.8840 1.76601 26
79         GP 7.9830 1.63145 27
80 Parametric 2.5022 0.31409 27
81       True 7.8840 1.76601 27
82         GP 6.3000 0.00000 28
83 Parametric 2.3670 0.23987 28
84       True 7.8840 1.76601 28
85         GP 7.9830 1.63145 29
86 Parametric 2.5945 0.23523 29
87       True 7.8840 1.76601 29
88         GP 8.0175 1.75271 30
89 Parametric 2.8178 0.08035 30
90       True 7.8840 1.76601 30
91         GP 8.1285 1.69568 31
92 Parametric 2.5992 0.27697 31
93       True 7.8840 1.76601 31
94         GP 6.3000 0.00000 32
95 Parametric 0.0660 0.11144 32
96       True 7.8840 1.76601 32
```

NULL