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

opts_knit$set(fig.path="figure/small-allee/")

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

load("small-allee.rda")

## ----figure_1, dependson=c("assemble-models", "par-fns"), 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"), 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."----
y <- numeric(8)
y[1] <- 4.5
for(t in 1:(length(y)-1))
      y[t+1] = z_g() * f(y[t], h=0, p=p)
# predicts means, does not reflect uncertainty estimate!
crash_data <- step_ahead_posteriors(y)
crash_data <- subset(crash_data, variable %in% c("GP", "Allen", "Ricker", "Myers"))
ggplot(crash_data) +
  geom_boxplot(aes(as.factor(as.integer(time)), value,
                   fill = variable, col=variable),
               alpha=.7, outlier.size=1, position="identity") +
#  geom_line(aes(time, value, col = variable,
#            group=interaction(L1,variable)), alpha=.1) +
  geom_point(aes(time, stock), size = 3) +
  scale_fill_manual(values=colorkey[c("GP", "Allen", "Ricker", "Myers")],
                      guide = guide_legend(override.aes = list(alpha = 1))) +
  scale_colour_manual(values=colorkey[c("GP", "Allen", "Ricker", "Myers")],
                      guide = guide_legend(override.aes = list(alpha = 1))) +
  facet_wrap(~variable) +
  theme(legend.position="none") + xlab("time") + ylab("stock size")

## ----figure_4, dependson = c("assemble-opt"), fig.cap="The steady-state optimal policy (infinite boundary) calculated under each model.  Policies are shown in terms of target escapement, $S_t$, as under models such as this a constant escapement policy is expected to be optimal [@Reed1979]."----
policies <- melt(data.frame(stock=x_grid, sapply(OPT, function(x) x_grid[x])), id="stock")
names(policies) <- c("stock", "method", "value")

ggplot(policies, aes(stock, stock - value, color=method)) +
  geom_line(lwd=1.2, alpha=0.8) + xlab("stock size") + ylab("escapement")  +
  scale_colour_manual(values=colorkey)

## ----figure_5, dependson=c("sim"), fig.cap="Gaussian process inference outperforms parametric estimates. Shown are 100 replicate simulations of the stock dynamics (eq 1) under the policies derived from each of the estimated models, as well as the policy based on the exact underlying model.", fig.width=6, fig.height=6----
ggplot(sims_data) +
  geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=.1) +
  scale_colour_manual(values=colorkey, guide = guide_legend(override.aes = list(alpha = 1))) +
  facet_wrap(~method) + guides(legend.position="none")

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

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

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

## ----figure_S1, fig.cap = "Traces from the MCMC estimates of the GP model show reasonable mixing (no trend) and sampling rejection rate (no piecewise jumps)"----
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)

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

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

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

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

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

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


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

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

```{r Ricker}
```

We use the Ricker model

```{r sdp-pars-ricker}
sigma_g <- 0.05
z_g <- function(sigma_g) rlnorm(1, 0, sigma_g) #1+(2*runif(1, 0,  1)-1)*sigma_g #
x_grid <- seq(0, 2.5 * K, length=101)
h_grid <- x_grid
profit = function(x,h) pmin(x, h)
delta <- 0.01
OptTime = 20
reward = profit(x_grid[length(x_grid)], x_grid[length(x_grid)]) + 1 / (1 - delta) ^ OptTime
xT <- 0
```

with parameters `r p`.

```{r x0-obs}
x_0_observed <- K
xT <- 0
set.seed(1)
```

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

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

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

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

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

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

Estimate data under the BH model

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

We fit a Gaussian process with

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

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

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

## The transition matrix of the inferred process

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

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

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

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

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

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

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

```{r policy_plot}
```

```{r othernoise}
```

```{r stationary_policy_only}
```

```{r simulate}
```

```{r tidy}
```

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

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

```{r costs}
```

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

repo_name: cboettig/nonparametric-bayes

repo_language: TeX

mime_type: text/plain


text:
`ro comment=NA, message=FALSE, cache.path="compare/", tidy=FALSE or`

## Sequential method

```{r}
require(MASS)
require(ggplot2)
require(kernlab)
opts_knit$set(upload.fun = socialR::notebook.url)
```

Parameterization-specific

```{r}
X <- seq(-5,5, length=40)
obs <- data.frame(x = c(-3, -1, 3),
                  y = c(0,  1, -1))
#obs <- data.frame(x = -5:5, y = sin(x) + rnorm(length(x),sd=.1))
l <- 1
sigma_n <- 0.5
```

Radial basis function/Gaussian kernel:

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

### Direct method

``` {r}
  n <- length(obs$x)
  K <- cov(obs$x, obs$x)
  I <-  diag(1, n)

  cov_xx_inv <- solve(K + sigma_n^2*I)
  Ef <- cov(X, obs$x) %*% cov_xx_inv %*% obs$y
  Cf <- cov(X, X) - cov(X, obs$x)  %*% cov_xx_inv %*% cov(obs$x, X)
```

### Direct sequential method, avoids matrix inverse instability

``` {r}

mmult <- function(x,y){
  if(length(x) == 1){
    x <- as.numeric(x)
    x * y
  } else if(length(y) == 1){
    y <- as.numeric(y)
    x * y
  }  else
  x %*% y
}

mu <- numeric(length(obs$y))
y <- obs$y
x <- obs$x

C_seq <- function(X, X_prime, i){
  if(i <= 1)
    cov(X, X_prime) - mmult(cov(X,x[i]), cov(x[i], X_prime)) / as.numeric( cov(x[i], x[i]) + sigma_n^2)
  else
    C_seq(X, X_prime,   i-1) - mmult(C_seq(X,x[i],   i-1), C_seq(x[i], X_prime,   i-1)) / as.numeric( C_seq(x[i], x[i],   i-1)  + sigma_n^2  )
}
mu_seq <- function(X, i){
  if(i <= 1)
    cov(x[i], X) * (y[i]-mu[i]) / as.numeric( cov(x[i], x[i]) + sigma_n^2)
  else
    mu_seq(X, i-1) + C_seq(x[i], X,  i-1) * (y[i]-mu[i]) / as.numeric( C_seq(x[i], x[i], i-1)  + sigma_n^2 )
}
ef <- t(mu_seq(X, length(obs$x)))
cf <- C_seq(X, X, length(obs$x))

```

Compare these results:

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


text:
```r
library(pomp)
```

```
## Loading required package: mvtnorm
```

```
## Loading required package: subplex
```

```
## Loading required package: deSolve
```

```r
opts_knit$set(upload.fun = socialR::flickr.url)
options(flickr_tags = "nonparametric-bayes")
```

```r
process.sim <- function(x, t, params, delta.t, ...) {
    # r = 0.75, K = 10, Q = 3, H = 1, a = 1.3
    r <- params["r"]
    K <- params["K"]
    Q <- params["Q"]
    H <- params["H"]
    a <- params["a"]
    sigma <- params["sigma"]
    X <- x["X"]
    epsilon <- exp(rlnorm(n = 1, mean = 0, sd = sigma))

    c(X = unname(X * exp(r * delta.t * (1 - X/K) - a * X^(Q - 1)/(X^Q + H^Q)) *
        epsilon))
}
```

```r
measure.sim <- function(x, t, params, ...) {
    tau <- params["tau"]
    X <- x["X"]
    y <- c(Y = unname(rlnorm(n = 1, meanlog = log(X), sd = tau)))
}
```

Create a container of class `pomp`

```r
may <- pomp(data = data.frame(time = 1:100, Y = NA), times = "time",
    rprocess = discrete.time.sim(step.fun = process.sim, delta.t = 1), rmeasure = measure.sim,
    t0 = 0)
```

Parameters and an initial condition

```r
theta <- c(r = 0.75, K = 10, Q = 3, H = 1, a = 1.3, sigma = 0.2,
    tau = 0.1, X.0 = 5)
```

Simulate

```r
dat <- simulate(may, params = theta, obs = TRUE, states = TRUE)
library(reshape2)
library(ggplot2)
dt <- melt(dat)[3:5]
names(dt) = c("time", "value", "variable")
ggplot(dt, aes(time, value, lty = variable)) + geom_line()
```

![plot of chunk unnamed-chunk-6](http://farm9.staticflickr.com/8527/8575983353_a8cbbdd4eb_o.png)

Note that simulate returns an instance of our object, with the data column filled in.  This is not true if we give it arguments such as state or obs = TRUE as above.  However, when we query for observed and true states in this way, the simulated data is not stored in our object.  We must assign it with `may <- pomp(data = ...` where data is a data frame with columns "time" and "Y".  This is much simpler if we use the default return:

```r
may <- simulate(may, params = theta)
```

Measurement density

```r
measure.density <- function(y, x, t, params, log, ...) {
    tau <- params["tau"]
    ## state at time t:
    X <- x["X"]
    ## observation at time t:
    Y <- y["Y"]
    ## compute the likelihood of Y|X,tau
    f <- dlnorm(x = Y, meanlog = log(X), sdlog = tau, log = log)

}
```

Add to our container

```r
may <- pomp(may, dmeasure = measure.density)
```

Optionally, work in the transformed variable space instead by providing parameter and inverse parameter transforms.

```r
may <- pomp(may, parameter.transform = function(params, ...) {
    exp(params)
}, parameter.inv.transform = function(params, ...) {
    log(params)
})
```

Run the particle filter.  Note that this evaluates the likelihood of our simulated data (stored in the object) at the given (true) parameters.

```r
pf <- pfilter(may, params = theta, Np = 1000)
logLik(pf)
```

```
## [1] -304.7
```

Here we go, estimate the model by particle filter.  This'll be slow...

```r
theta.guess <- theta + c(abs(rnorm(length(theta)-1)),0)
estpars <- names(theta[1:7]) # which pars should be estimated. lets try and get them all
rw.sd <- theta[1:7]/(50*theta[1:7]) # a named vector of random walk values. This sets it to 0.02 for all pars
system.time(mf <-
  mif(may, Nmif = 10, # Number of iterations
      start = theta.guess, # Starting value of parameters
      transform = TRUE, # work in transformed parameter space (e.g. logs avoid neg par values)
      pars = estpars,  # which parameters to estimate (others fixed at starting value)
      rw.sd = rw.sd, # random walk width
      Np = 200, # Number of particles
      var.factor = 4, # scaling factor for random walk sd
      ic.lag = 10, # fixed lag
      cooling.factor = 0.999, # exponential cooling
      max.fail = 10) # internal tolerance
)
```

```
##    user  system elapsed
##  22.018   0.008  22.106
```

```r
compare.mif(mf)
```

![plot of chunk unnamed-chunk-13](http://farm9.staticflickr.com/8507/8575984095_b391cbb8c8_o.png) ![plot of chunk unnamed-chunk-13](http://farm9.staticflickr.com/8382/8577082802_82dc557940_o.png)

```r
coef(mf)
```

```
##       r       K       Q       H       a   sigma     tau     X.0
##  0.6985 10.2655  4.0470  1.2844  1.6782  0.1587  0.1711  5.0000
```

```r
logLik(mf)
```

```
## [1] -314.4
```

```r
pf_est <- pfilter(may, params = theta.guess, Np = 1000)
logLik(pf_est)  # initial loglik
```

```
## [1] -441.9
```

These functions apply to replicated mifs above.

```r
theta.true <- coef(may)
theta.mif <- apply(sapply(mf, coef), 1, mean)
```

```
## Error: no method for coercing this S4 class to a vector
```

```r
loglik.mif <- replicate(n = 10, logLik(pfilter(mf[[1]], params = theta.mif,
    Np = 10000)))
```

```
## recover called non-interactively; frames dumped, use debugger() to view
## recover called non-interactively; frames dumped, use debugger() to view
```

```
## Error: error in evaluating the argument 'object' in selecting a method for
## function 'logLik': Error in pfilter(mf[[1]], params = theta.mif, Np =
## 10000) : error in evaluating the argument 'object' in selecting a method
## for function 'pfilter': Error in mf[[1]] : this S4 class is not
## subsettable
```

```r
bl <- mean(loglik.mif)
```

```
## Error: object 'loglik.mif' not found
```

```r
loglik.mif.est <- bl + log(mean(exp(loglik.mif - bl)))
```

```
## Error: object 'bl' not found
```

```r
loglik.mif.se <- sd(exp(loglik.mif - bl))/exp(loglik.mif.est - bl)
```

```
## Error: object 'loglik.mif' not found
```


text:
```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
```

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

![plot of chunk unnamed-chunk-7](http://carlboettiger.info/assets/figures/2012-12-07-c431834819-unnamed-chunk-7.png)

### kernlab method

```r
k <- gp_fit(obs, X, c(sigma_n=sigma_n, l=l, tau=1), "kernlab")
```

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

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

![plot of chunk unnamed-chunk-9](http://carlboettiger.info/assets/figures/2012-12-07-c431834819-unnamed-chunk-9.png)


text:
```{r}
library(pomp)
#opts_knit$set(upload.fun = socialR::flickr.url)
#options(flickr_tags = "nonparametric-bayes")
```

```{r}
process.sim <- function(x, t, params, delta.t, ...){
# r = 0.75, K = 10, Q = 3, H = 1, a = 1.3
  r <- params["r"]
  K <- params["K"]
  Q <- params["Q"]
  H <- params["H"]
  a <- params["a"]
  sigma <- params["sigma"]
  X <- x["X"]
  epsilon <- exp(rlnorm(n = 1, mean = 0, sd = sigma))

  c(X = unname(
    X * exp(r * delta.t * (1 - X / K) - a * X ^ (Q - 1) / (X ^ Q + H ^ Q)) * epsilon))
}
```

```{r}
measure.sim <- function(x, t, params, ...){
  tau <- params["tau"]
  X <- x["X"]
  y <- c(Y = unname(rlnorm(n = 1, meanlog = log(X), sd=tau)))
}
```

Create a container of class `pomp`

```{r}
may <- pomp(data = data.frame(time = 1:100, Y = NA), times = "time", rprocess = discrete.time.sim(step.fun = process.sim, delta.t = 1), rmeasure = measure.sim, t0 = 0)
```

Parameters and an initial condition

```{r}
theta <- c(r = 0.75, K = 10, Q = 3, H = 1, a = 1.3,
           sigma = 0.2, tau = 0.1, X.0 = 5)
```

Simulate

```{r}
dat <- simulate(may, params=theta, obs=TRUE, states=TRUE)
library(reshape2)
library(ggplot2)
dt <- melt(dat)[3:5]
names(dt) = c("time", "value", "variable")
ggplot(dt, aes(time, value, lty=variable)) + geom_line()
```

Note that simulate returns an instance of our object, with the data column filled in.  This is not true if we give it arguments such as state or obs = TRUE as above.  However, when we query for observed and true states in this way, the simulated data is not stored in our object.  We must assign it with `may <- pomp(data = ...` where data is a data frame with columns "time" and "Y".  This is much simpler if we use the default return:

```{r}
may <- simulate(may, params=theta)
```

Measurement density

```{r}
measure.density <- function(y, x, t, params, log, ...) {
    tau <- params["tau"]
    ## state at time t:
    X <- x["X"]
    ## observation at time t:
    Y <- y["Y"]
    ## compute the likelihood of Y|X,tau
    f <- dlnorm(x = Y, meanlog = log(X), sdlog = tau, log = log)

}
```

Add to our container

```{r}
may <- pomp(may, dmeasure = measure.density)
```

Optionally, work in the transformed variable space instead by providing parameter and inverse parameter transforms.

```{r}
may  <- pomp(may, parameter.transform = function(params,
    ...) {
    exp(params)
}, parameter.inv.transform = function(params, ...) {
    log(params)
})
```

Run the particle filter.  Note that this evaluates the likelihood of our simulated data (stored in the object) at the given (true) parameters.

```{r}
pf <- pfilter(may, params = theta, Np = 1000)
logLik(pf)
```

Here we go, estimate the model by particle filter.  This'll be slow...

```{r}
theta.guess <- theta + c(abs(rnorm(length(theta)-1)),0)
estpars <- names(theta[1:7]) # which pars should be estimated. lets try and get them all
rw.sd <- theta[1:7]/(50*theta[1:7]) # a named vector of random walk values. This sets it to 0.02 for all pars
system.time(mf <-
  mif(may, Nmif = 10, # Number of iterations
      start = theta.guess, # Starting value of parameters
      transform = TRUE, # work in transformed parameter space (e.g. logs avoid neg par values)
      pars = estpars,  # which parameters to estimate (others fixed at starting value)
      rw.sd = rw.sd, # random walk width
      Np = 200, # Number of particles
      var.factor = 4, # scaling factor for random walk sd
      ic.lag = 10, # fixed lag
      cooling.factor = 0.999, # exponential cooling
      max.fail = 10) # internal tolerance
)
```

```{r}
compare.mif(mf)
```

```{r}
coef(mf)
logLik(mf)
pf_est <- pfilter(may, params = theta.guess, Np = 1000)
logLik(pf_est) # initial loglik
```

These functions apply to replicated mifs above.

```{r}
theta.true <- coef(may)
theta.mif <- apply(sapply(mf,coef),1,mean)
loglik.mif <- replicate(n=10,logLik(pfilter(mf[[1]],params=theta.mif,Np=10000)))
bl <- mean(loglik.mif)
loglik.mif.est <- bl+log(mean(exp(loglik.mif-bl)))
loglik.mif.se <- sd(exp(loglik.mif-bl))/exp(loglik.mif.est-bl)
```


text:
multi-dimensional dynamics example
========================================================

```r
library(reshape2)
library(ggplot2)
library(data.table)
library(kernlab)
library(tgp)
```

$$\begin{align}
\frac{dM}{dt} = aMC - \frac{(g-h) M}{M+T} + \gamma M T \\
\frac{dC}{dt} = rTC - dC - a M C
\end{align}$$

```r
mumby <- function(x, h, p= c(a = 0.1, g = .5, gamma = 0.8, r = 1, d = 0.44)){
  M <- x[1] # Macroalgae
  C <- x[2] # Corals
  a <- p[1] # algal growth rate
  g <- p[2] # grazing rate (reduced by harvesting)
  gamma <- p[3] # algal colinization of dead tufts
  r <- p[4] # Coral growth rate
  d <- p[5] # Coral death rate
  dt <- 0.025 #
  M_t <- M * exp( dt * (a * M * C - (g - h) * M / (M + (1-M-C)) + gamma * M * (1-M-C)) )
  C_t <- C * exp( dt * (r * (1-M-C) * C - d * C - a * M * C) )
  c(M_t, C_t)
}

```

```r
Tmax <- 2500
X <- matrix(NA, nrow=Tmax, ncol=2)
colnames(X) = c("M", "C")
```

```r
X[1,] = c(0.9,0.9)
for(i in 1:(Tmax-1))
  X[i+1,] = mumby(X[i,],0)
dat <- data.frame(time = 1:Tmax, X)
dat <- dat[seq(1,2500,by=25),]
```

```r
ggplot(dat) + geom_line(aes(time, M), col="red") + geom_line(aes(time, C), col="blue")
```

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

```r
inits <- as.matrix(expand.grid(seq(0.05,0.6,length=11), seq(0.05,0.6, length=11)))
df <- lapply(1:round(length(inits[,1])), function(j){
  X[1,] = inits[j,]
  for(i in 1:(Tmax-1))
    X[i+1,] = mumby(X[i,],0.2)
  data.frame(time = 1:Tmax, X)
  })
df2 <- melt(df, id=c("time", "M", "C"))
ggplot(df2) + geom_line(aes(M, C, group=L1)) +
  geom_point(data = subset(df2, time==1), aes(M, C, group=L1), col="red")
```

![plot of chunk phase-space](http://carlboettiger.info/assets/figures/2012-12-11-b74ccc4566-phase-space.png)

Noisy data example

```r
sigma_g = 0.1
X[1,] = c(0.9,0.9)
for(i in 1:(Tmax-1))
  X[i+1,] = rlnorm(1,0, sigma_g) * mumby(X[i,],0)
dat <- data.frame(time = 1:Tmax, X)
dat <- dat[seq(1,2500,by=25),]
```

Multi-dimensional GP, `tgp` method, based on time-delay data:

```r
n <- length(dat[,1])
obs <- data.frame(x1=dat[2:(n-1),2], x2 = dat[1:(n-2),2], y=dat[3:n,2])
x_grid <- seq(0,1,length.out=100)
```

```r
require(rgl)
```

```
Loading required package: rgl
```

```r
plot3d(dat$time, dat$M, dat$C)
plot3d(obs$x1, obs$x2, obs$y)
```

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

```r
plot(gp, layout="surf")
```

![plot of chunk unnamed-chunk-8](http://carlboettiger.info/assets/figures/2012-12-11-b74ccc4566-unnamed-chunk-8.png)

Extract the posterior Gaussian process mean and the $\pm 2$ standard deviations over the predicted grid from the fit:

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

```r
ggplot(tgp_dat)  + geom_ribbon(aes(x1,y,ymin=ymin,ymax=ymax), fill="gray80") +
  geom_line(aes(x1,y)) + geom_point(data=obs, aes(x1,y))
```

![plot of chunk x1](http://carlboettiger.info/assets/figures/2012-12-11-b74ccc4566-x1.png)

```r
ggplot(tgp_dat)  + geom_ribbon(aes(x2,y,ymin=ymin,ymax=ymax), fill="gray80") +
  geom_line(aes(x2,y)) + geom_point(data=obs, aes(x2,y))
```

![plot of chunk x2](http://carlboettiger.info/assets/figures/2012-12-11-b74ccc4566-x2.png)

`kernlab` method


text:
```r
library(MCMCpack)
library(tgp)
library(reshape2)
library(ggplot2)
```

I am trying to understand the interface for tgp method, but much of it is not well documented (yes, despite reading through the package R manual, two nice vignettes, and a the PhD Thesis describing it.)

Just to get some consistent notation down: consider the case of a Gaussian process with a Gaussian/radial basis function kernel, parameterized as:

$$K(X, X') = \sigma^2 e^{\frac{(X-X')^2}{d}}$$

And observations from $Z = f(X) + \varepsilon$, for which we seek to approximate $f$ as a (multivariate) Gaussian process, $Z | X ~ N(\mu, C)$. The covariance matrix $C$ is given by our kernel $K$, conditioned on our observations; the mean $\mu$ is given by a constant or linear model, likewise conditioned on our observations.  If the multivariate normal of observed and predicted points 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)$$

the conditional probability is

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

Our hyperparameters are $d$, $\sigma^2$, and $\varepsilon$, for which we need priors.  We'll use the inverse gamma,

$$P(x; a, g) = \frac{g^a}{\Gamma(a)} x^{-a - 1}\exp\left(-\frac{g}{x}\right)$$

The distribution arises as the marginal posterior distribution for the unknown variance of a normal distribution if an uninformative prior is used; and as a useful conjugate prior if a less uninformative prior is preferred. However, it is common among Bayesians to consider an alternative parametrization of the normal distribution in terms of the precision, defined as the reciprocal of the variance, which allows the gamma distribution to be used directly as a conjugate prior.

* $X$ is `X`, $Z$ is `Z`

* $\sigma^2$ is the parameter `s2` in the `tgp` package, coming from inverse gamma prior with parameters $a = a_0$ and $g = g_0$, which are set by `s2.p = c(a_0, g_0)`.  Optionally, $a_0$ and $g_0$ can be taken from an exponential prior $\lambda e^{-\lambda x}$, in which the value of $\lambda$ is given by `s2.lam = c(lambda_a, lambda_g)`.  Setting `s2.lam = "fixed"` "turns off" the use of a hyperprior.

* $d$ I believe is `d`, which Grammarcy calls the range parameter (other terms include the "lengthscale"). I am not sure what either `d` or `s2` correspond to if the covariance function `corr` is something other than the radial basis function.  I'm also not clear how to set the correlation function to a strictly Gaussian covariance -- `corr` is set to the `exp` family, which leaves the power as a parameter `p_0`:

$$K(X, X') =  e^{\frac{-(X-X')^p_0}{d}}$$

(from equation 4 in the Vignette (1)).  The parameter `p_0` is fixed at `2`, and not treated as a hyperparameter.

The length-scale comes from a hyperparameter that is the sum (mixture) of two Gamma distributions, with prior parameter `c(a1,g1,a2,g2)` where `g1` and `g2` are scale (1/rate) parameters.  As before, these parameters  `c(a1, g1, a2, g2)` can be drawn from an exponential distribution hyperprior whose $\lambda$ values are given by `d.lam`, or turned off.

* $\varepsilon$ is `nug`.  It seems to have same shape prior as $d$, specified by `c(a1,g1,a2,g2)`, and optional hyperprior as well.  Additionally it can be fixed to a constant by setting `nug.p = 0` and specifying the staring value as `gd[1]`.

* The mean is not fixed at zero but is instead a linear model, we need hyperparameters for it as well.  The key parameter there is the slope $\beta$.  The functional form for it's prior is specified as `bprior`, which takes values `bflat`, `b0`, `bmzt`, or `bmznot`.  A sensible choice is as a Gaussian `b0`, as shown in Equation 1 of the vignette

$$
\begin{align}
Z | \beta, \sigma^2, K &\sim N_n(\mathbf{F} \beta, \sigma^2 \mathbf{K}) \\
\beta | \sigma^2, \tau^2, \mathbf{W}, \beta_0 &\sim N_m(\beta_0, \sigma^2 \tau^2 \mathbf{W} )  \\
\beta_0 &\sim N_m(\mathbf{\mu}, \mathbf{B})
\end{align}
$$

## Example comparison

Say, given some observed data pairs $x,y$ and predictions desired on a grid $X$,

```r
X = seq(-5, 5, length=50)
x = c(-4, -3, -1,  0,  2)
y = c(-2,  0,  1,  2, -1)
```

Manually I could just do:

```r
d = .5; epsilon = .1; sigma = 1; #fixed hyperparamaters
SE <- function(Xi,Xj, d) sigma * exp(-0.5 * (Xi - Xj) ^ 2 / d ^ 2)
cov <- function(X, Y) outer(X, Y, SE, d)
cov_xx_inv <- solve(cov(x, x) + epsilon * diag(1, length(x)))
Ef <- cov(X, x) %*% cov_xx_inv %*% y
Cf <- cov(X, X) - cov(X, x)  %*% cov_xx_inv %*% cov(x, X)
```

And `Ef` and `Cf` are the estimated (krieging) mean and covariance.

## Priors

Note: it appears that the mixed gammas for `d` and `nug` actually treat the second argument, `g1` and `g2`, as the rate instead of the scale, so I invert the desired scale:

```r
s2.p <- c(50,50)
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)
```

Define as curves and plot the priors:

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

xx <- seq(.0001, 5, length.out=100)
priors <- data.frame(x = xx, s2 = s2_prior(xx), tau2 = tau2_prior(xx), beta0 = beta0_prior(xx, 1), nug = nug_prior(xx), d = d_prior(xx))
priors <- melt(priors, id="x")
ggplot(priors) + geom_line(aes(x, value)) + facet_wrap(~variable, scale="free")
```

![plot of chunk unnamed-chunk-5](http://carlboettiger.info/assets/figures/2012-12-10-9b2d792442-unnamed-chunk-5.png)

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

Extract the posterior Gaussian process mean and the $\pm 2$ standard deviations over the predicted grid from the fit:

```r
V <- gp$ZZ.ks2
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))
```

and likewise for our manual data:

```r
manual_dat <- data.frame(x = X,
                         y = Ef,
                         ymin = (Ef - 1.96 * sqrt(diag(Cf))),
                         ymax = (Ef + 1.96 * sqrt(diag(Cf))))
```

Compare the GP posteriors:

```r
ggplot(tgp_dat) +
    geom_ribbon(aes(x, y, ymin = ymin, ymax = ymax), fill="red", alpha = .1) + # Var
    geom_line(aes(x, y), col="red") + # mean
    geom_point(data = data.frame(x = x, y = y), aes(x, y)) + # raw data
    geom_ribbon(dat = manual_dat, aes(x, y, ymin = ymin, ymax = ymax), fill = "blue", alpha = .1) + # Var
    geom_line(dat = manual_dat, aes(x, y), col = "blue")  + #MEAN
    theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA))
```

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

## Posteriors of hyperparameters

```r
require(reshape2)
hyperparameters <- c("index", "s2", "tau2", "beta0", "nug", "d", "ldetK")
posteriors <- melt(gp$trace$XX[[1]][,hyperparameters], id="index")
ggplot(posteriors) + geom_histogram(aes(value)) + facet_wrap(~variable, scales="free")
```

![plot of chunk unnamed-chunk-9](http://carlboettiger.info/assets/figures/2012-12-10-9b2d792442-unnamed-chunk-9.png)


text:
Exploring the behavior and ifluence of the prior in the `tgp` package.  It appears there may be a bug (or rather incorrect documentation) in how the mixed Gamma priors are defined for the nugget `nug` and range / length-scale `d`.

```{r include=FALSE}
opts_knit$set(upload.fun = socialR::flickr.url)
opts_chunk$set(dev.args=list(bg="transparent"), comment=NA, tidy=FALSE, message=FALSE)
```

```{r libraries}
library(MCMCpack)
library(tgp)
library(reshape2)
library(ggplot2)
```

Consider some very simple data:

```{r}
X = seq(-5, 5, length=50)
x = c(-4, -3, -1,  0,  2)
y = c(-2,  0,  1,  2, -1)
```

## Priors

We choose prior parameters that give us very strong priors around particular values.

```{r}
s2.p <- c(50,50) # InvGamma, approaches Gaussian approaching delta-fn around 1
tau2.p <- c(20,1) # InvGamma, approaches exponential approaching delta-fn around 0
d.p = c(10, 0.01, 10, 0.01) # Mixed Gammas, approaches Gaussian around 0.1, with var .001
nug.p = c(10, 0.1, 10, 0.1) # Mixed Gammas, approaches Gaussian around 1, with var .1
```

Visualize the priors on lengthscale and nugget:

```{r}
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])
xx <- seq(.0001, 2, length.out=100)
priors <- data.frame(x = xx, nug = nug_prior(xx), d = d_prior(xx))
priors <- melt(priors, id="x")
ggplot(priors) + geom_line(aes(x, value)) + facet_wrap(~variable, scale="free")
```

Solve the gp with these priors:

```{r}
gp <- bgp(X=x, XX=X, Z=y, verb=0,
          meanfn="constant", bprior="b0", BTE=c(1,2000,1), m0r1=FALSE,
          corr="exp", trace=TRUE, beta = 0,
          s2.p = s2.p, d.p = d.p, nug.p = nug.p,
          s2.lam = "fixed", d.lam = "fixed", nug.lam = "fixed",
          tau2.lam = "fixed", tau2.p = tau2.p)
```

Extract the posterior Gaussian process mean and the $\pm 2$ standard deviations over the predicted grid from the fit:

```{r}
V <- gp$ZZ.ks2
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))
```

For comparison, let us manually estimate the GP with hyperparameters fixed at the mean value of the priors.

```{r}
d = .1; epsilon = .1; sigma = 1; #fixed hyperparamaters
SE <- function(Xi,Xj, d) sigma * exp(- (Xi - Xj) ^ 2 / d)
cov <- function(X, Y) outer(X, Y, SE, d)
cov_xx_inv <- solve(cov(x, x) + epsilon * diag(1, length(x)))
Ef <- cov(X, x) %*% cov_xx_inv %*% y
Cf <- cov(X, X) - cov(X, x)  %*% cov_xx_inv %*% cov(x, X)
manual_dat <- data.frame(x = X,
                         y = Ef,
                         ymin = (Ef - 1.96 * sqrt(diag(Cf))),
                         ymax = (Ef + 1.96 * sqrt(diag(Cf))))
```

Compare the GP posteriors in a plot:

```{r gp-plot}
ggplot(tgp_dat) +
    geom_ribbon(aes(x, y, ymin = ymin, ymax = ymax), fill="red", alpha = .1) + # Var
    geom_line(aes(x, y), col="red") + # mean
    geom_point(data = data.frame(x = x, y = y), aes(x, y)) + # raw data
    geom_ribbon(dat = manual_dat, aes(x, y, ymin = ymin, ymax = ymax), fill = "blue", alpha = .1) + # Var
    geom_line(dat = manual_dat, aes(x, y), col = "blue")  + #MEAN
    theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA))
```

Wait a minute. The Bayesian estimate of the lengthscale seems way too long, given our strongly informative prior. Likewise the measurement uncertainty appears to have induced a rather broad distribution.  Let's look at the posteriors:

```{r}
require(reshape2)
hyperparameters <- c("index", "s2", "tau2", "beta0", "nug", "d", "ldetK")
posteriors <- melt(gp$trace$XX[[1]][,hyperparameters], id="index")
ggplot(posteriors) + geom_histogram(aes(value)) + facet_wrap(~variable, scales="free")

```

Hmm... posteriors indeed show a huge value for the length scale `d`, and for the nugget `nug`.  Looks supiciously like the mixed Gamma prior when the scale parameter is used as the rate parameter (it's inverse):

## Compare to priors with `scale` misplaced as `rate` in `d`, `nug`

```{r density_prior}
d_prior <- function(x) dgamma(x, d.p[1], d.p[2]) + dgamma(x, d.p[3], d.p[4])
nug_prior <- function(x) dgamma(x, nug.p[1], nug.p[2]) + dgamma(x, nug.p[3], nug.p[4])
xx <- seq(.0001, 2500, length.out=1000)
priors <- data.frame(x = xx, nug = nug_prior(xx), d = d_prior(xx))
priors <- melt(priors, id="x")
ggplot(priors) + geom_line(aes(x, value)) + facet_wrap(~variable, scale="free")
```


text:
Exploring the behavior and ifluence of the prior in the `tgp` package.  It appears there may be a bug (or rather incorrect documentation) in how the mixed Gamma priors are defined for the nugget `nug` and range / length-scale `d`.

```r
library(MCMCpack)
library(tgp)
library(reshape2)
library(ggplot2)
```

Consider some very simple data:

```r
X = seq(-5, 5, length=50)
x = c(-4, -3, -1,  0,  2)
y = c(-2,  0,  1,  2, -1)
```

## Priors

We choose prior parameters that give us very strong priors around particular values.

```r
s2.p <- c(50,50) # InvGamma, approaches Gaussian approaching delta-fn around 1
tau2.p <- c(20,1) # InvGamma, approaches exponential approaching delta-fn around 0
d.p = c(10, 0.01, 10, 0.01) # Mixed Gammas, approaches Gaussian around 0.1, with var .001
nug.p = c(10, 0.1, 10, 0.1) # Mixed Gammas, approaches Gaussian around 1, with var .1
```

Visualize the priors on lengthscale and nugget:

```r
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])
xx <- seq(.0001, 2, length.out=100)
priors <- data.frame(x = xx, nug = nug_prior(xx), d = d_prior(xx))
priors <- melt(priors, id="x")
ggplot(priors) + geom_line(aes(x, value)) + facet_wrap(~variable, scale="free")
```

![plot of chunk unnamed-chunk-4](http://farm9.staticflickr.com/8404/8694391470_8ac745ba0f_o.png)

Solve the gp with these priors:

```r
gp <- bgp(X=x, XX=X, Z=y, verb=0,
          meanfn="constant", bprior="b0", BTE=c(1,2000,1), m0r1=FALSE,
          corr="exp", trace=TRUE, beta = 0,
          s2.p = s2.p, d.p = d.p, nug.p = nug.p,
          s2.lam = "fixed", d.lam = "fixed", nug.lam = "fixed",
          tau2.lam = "fixed", tau2.p = tau2.p)
```

Extract the posterior Gaussian process mean and the $\pm 2$ standard deviations over the predicted grid from the fit:

```r
V <- gp$ZZ.ks2
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))
```

For comparison, let us manually estimate the GP with hyperparameters fixed at the mean value of the priors.

```r
d = .1; epsilon = .1; sigma = 1; #fixed hyperparamaters
SE <- function(Xi,Xj, d) sigma * exp(- (Xi - Xj) ^ 2 / d)
cov <- function(X, Y) outer(X, Y, SE, d)
cov_xx_inv <- solve(cov(x, x) + epsilon * diag(1, length(x)))
Ef <- cov(X, x) %*% cov_xx_inv %*% y
Cf <- cov(X, X) - cov(X, x)  %*% cov_xx_inv %*% cov(x, X)
manual_dat <- data.frame(x = X,
                         y = Ef,
                         ymin = (Ef - 1.96 * sqrt(diag(Cf))),
                         ymax = (Ef + 1.96 * sqrt(diag(Cf))))
```

Compare the GP posteriors in a plot:

```r
ggplot(tgp_dat) +
    geom_ribbon(aes(x, y, ymin = ymin, ymax = ymax), fill="red", alpha = .1) + # Var
    geom_line(aes(x, y), col="red") + # mean
    geom_point(data = data.frame(x = x, y = y), aes(x, y)) + # raw data
    geom_ribbon(dat = manual_dat, aes(x, y, ymin = ymin, ymax = ymax), fill = "blue", alpha = .1) + # Var
    geom_line(dat = manual_dat, aes(x, y), col = "blue")  + #MEAN
    theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA))
```

![plot of chunk gp-plot](http://farm9.staticflickr.com/8405/8693273719_aa17e176c7_o.png)

Wait a minute. The Bayesian estimate of the lengthscale seems way too long, given our strongly informative prior. Likewise the measurement uncertainty appears to have induced a rather broad distribution.  Let's look at the posteriors:

```r
require(reshape2)
hyperparameters <- c("index", "s2", "tau2", "beta0", "nug", "d", "ldetK")
posteriors <- melt(gp$trace$XX[[1]][,hyperparameters], id="index")
ggplot(posteriors) + geom_histogram(aes(value)) + facet_wrap(~variable, scales="free")
```

![plot of chunk unnamed-chunk-8](http://farm9.staticflickr.com/8115/8693274145_a84a10948c_o.png)

Hmm... posteriors indeed show a huge value for the length scale `d`, and for the nugget `nug`.  Looks supiciously like the mixed Gamma prior when the scale parameter is used as the rate parameter (it's inverse):

## Compare to priors with `scale` misplaced as `rate` in `d`, `nug`

```r
d_prior <- function(x) dgamma(x, d.p[1], d.p[2]) + dgamma(x, d.p[3], d.p[4])
nug_prior <- function(x) dgamma(x, nug.p[1], nug.p[2]) + dgamma(x, nug.p[3], nug.p[4])
xx <- seq(.0001, 2500, length.out=1000)
priors <- data.frame(x = xx, nug = nug_prior(xx), d = d_prior(xx))
priors <- melt(priors, id="x")
ggplot(priors) + geom_line(aes(x, value)) + facet_wrap(~variable, scale="free")
```

![plot of chunk density_prior](http://farm9.staticflickr.com/8393/8693274295_a568f0da02_o.png)


text:
```r
f <- Ricker
p <- c(2.5, 10)
K <- 12.5
```

We use the Ricker model

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

with parameters `2.5, 10`.

```r
x_0_observed <- K
xT <- 0
set.seed(1)
```

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

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

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

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

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

```r
estf <- function(p){
  mu <- 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']
```

```r
s2.p <- c(50,50)
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)
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)
```

Estimate data under the BH model

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

```
Warning: for memory/storage reasons, trace not recommended when
3*(10+d)*(BTE[2]-BTE[1])*R*(nn+1)/BTE[3]=23562000 > 1e+7.  Try reducing
nrow(XX)
```

We fit a Gaussian process with

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

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

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

```r
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 = dnorm, 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")
```

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

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

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

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

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

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

![plot of chunk gp-posteriors](http://carlboettiger.info/assets/figures/2012-12-19-4ce05b9191-gp-posteriors1.png)

```r
ggplot(prior_curves) +
  geom_line(aes(x=value, y=density), col="red", lwd=2) +
  facet_wrap(~ variable, scale="free")
```

![plot of chunk gp-posteriors](http://carlboettiger.info/assets/figures/2012-12-19-4ce05b9191-gp-posteriors2.png)

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

## The transition matrix of the inferred process

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

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

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

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

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

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

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

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

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

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

```r
policies <- melt(data.frame(stock=x_grid,
                            GP = x_grid[opt_gp$D[,1]],
                            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
```

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

```r
z_g = function() rlnorm(1, 0, sigma_g)
z_m = function() 1+(2*runif(1, 0,  1)-1) * sigma_m
```

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

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

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

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

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

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

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

```r
profits <- dt[, sum(profit), by = c("reps", "method")]
means <- profits[, mean(V1), by = method]
sds <- profits[, sd(V1), by = method]
cbind(means, sd = sds$V1)
```

```
       method     V1    sd
1:         GP  25.29 5.683
2: Parametric  12.50 0.000
3:       True 266.96 3.542
```

Myers RA, Barrowman NJ, Hutchings JA and Rosenberg AA (1995).
“Population Dynamics of Exploited Fish Stocks at Low Population Levels.”
Science, 269.
ISSN 0036-8075, http://dx.doi.org/10.1126/science.269.5227.1106.


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

```r
require(pdgControl)
require(nonparametricbayes)
require(ggplot2)
```

### Beverton-Holt function

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

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

Parameter definitions

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

Noise function, profit function

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

Simulation

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

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

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


text:
# Semi-analytic posteriors

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

### Generating Model and parameters

Ricker model, parameterized as

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

for random unit normal $Z_t$

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

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

### Sample Data

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

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

$$P(\beta | X) $$

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

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

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

beta_grid[which.min(P_B.X)]

```

# Estimating the Allen model:

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

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

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

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

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

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

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

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

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

and loop over the grid on each

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

````

We summarize with a contour plot

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

```

```{r include=FALSE}
# Conditional Distributions
tmpA <- sapply(a_grid, prob, b=0.1)
tmpB <- sapply(b_grid, function(b) prob(a=0.01, b=b))

# Modes
a_grid[which.min(tmpA)]
b_grid[which.min(tmpB)]

qplot(a_grid, -log(tmpA))
qplot(b_grid, -log(tmpB))
```

# Estimating the Myers model on this data:

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

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

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

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

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

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

```

```{r include=FALSE}
# Conditional Distributions
tmpK <- sapply(K_grid, function(k) prob(1,k))
tmptheta <- sapply(theta_grid, prob, K=10)

# Modes
K_grid[which.min(tmpK)]
theta_grid[which.min(tmptheta)]

# Means
K_grid %*% -log(tmpK) / sum( -log(tmpK) )
theta_grid %*% -log(tmptheta) / sum( -log(tmptheta) )

qplot(theta_grid, -log(tmptheta))
qplot(K_grid, -log(tmpK))
```

## Reconstructing the joint probability distribution

$$P(r, \beta, \sigma | X) = \frac{P(X | r, \beta, \sigma) P(r) P(\beta) P(\sigma) }{\int_r \int_{\beta} \int_{\sigma P(X | r, \beta, \sigma) P(r) P(\beta) P(\sigma)}  $$

Having done the first two integrals analytically and then evaluating the remaining expression over the grid, we simply sum over the grids to complete our calculation of the normalization term from the denominator,

```{r}
N_ricker <- sum(P_B.X) * (beta_grid[2]-beta_grid[1])
```

We then have an analytic function for the numerator,

$$ P(X | r, \beta, \sigma) P(r) P(\beta) P(\sigma) =  \frac{1}{\sqrt{2 \pi \sigma^2}^{T-1}} \exp\left(\frac{\sum_t^{T-1} \left(\log X_{t+1} - \log X_t - \log r + \beta X_t\right)^2 }{2 \sigma^2}\right) \frac{1}{\sigma^2}$$

```{r}
ricker_posterior <- function(r, beta, sigma){}
```

## Transforming back


text:
```{r include=FALSE}
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 libraries}
library(MCMCpack)
library(tgp)
library(reshape2)
library(ggplot2)
```

I am trying to understand the interface for tgp method, but much of it is not well documented (yes, despite reading through the package R manual, two nice vignettes, and a the PhD Thesis describing it.)

Just to get some consistent notation down: consider the case of a Gaussian process with a Gaussian/radial basis function kernel, parameterized as:

$$K(X, X') = \sigma^2 e^{\frac{(X-X')^2}{d}}$$

And observations from $Z = f(X) + \varepsilon$, for which we seek to approximate $f$ as a (multivariate) Gaussian process, $Z | X ~ N(\mu, C)$. The covariance matrix $C$ is given by our kernel $K$, conditioned on our observations; the mean $\mu$ is given by a constant or linear model, likewise conditioned on our observations.  If the multivariate normal of observed and predicted points 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)$$

the conditional probability is

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

Our hyperparameters are $d$, $\sigma^2$, and $\varepsilon$, for which we need priors.  We'll use the inverse gamma,

$$P(x; a, g) = \frac{g^a}{\Gamma(a)} x^{-a - 1}\exp\left(-\frac{g}{x}\right)$$

The distribution arises as the marginal posterior distribution for the unknown variance of a normal distribution if an uninformative prior is used; and as a useful conjugate prior if a less uninformative prior is preferred. However, it is common among Bayesians to consider an alternative parametrization of the normal distribution in terms of the precision, defined as the reciprocal of the variance, which allows the gamma distribution to be used directly as a conjugate prior.

* $X$ is `X`, $Z$ is `Z`

* $\sigma^2$ is the parameter `s2` in the `tgp` package, coming from inverse gamma prior with parameters $a = a_0$ and $g = g_0$, which are set by `s2.p = c(a_0, g_0)`.  Optionally, $a_0$ and $g_0$ can be taken from an exponential prior $\lambda e^{-\lambda x}$, in which the value of $\lambda$ is given by `s2.lam = c(lambda_a, lambda_g)`.  Setting `s2.lam = "fixed"` "turns off" the use of a hyperprior.

* $d$ I believe is `d`, which Grammarcy calls the range parameter (other terms include the "lengthscale"). I am not sure what either `d` or `s2` correspond to if the covariance function `corr` is something other than the radial basis function.  I'm also not clear how to set the correlation function to a strictly Gaussian covariance -- `corr` is set to the `exp` family, which leaves the power as a parameter `p_0`:

$$K(X, X') =  e^{\frac{-(X-X')^p_0}{d}}$$

(from equation 4 in the Vignette (1)).  The parameter `p_0` is fixed at `2`, and not treated as a hyperparameter.

The length-scale comes from a hyperparameter that is the sum (mixture) of two Gamma distributions, with prior parameter `c(a1,g1,a2,g2)` where `g1` and `g2` are scale (1/rate) parameters.  As before, these parameters  `c(a1, g1, a2, g2)` can be drawn from an exponential distribution hyperprior whose $\lambda$ values are given by `d.lam`, or turned off.

* $\varepsilon$ is `nug`.  It seems to have same shape prior as $d$, specified by `c(a1,g1,a2,g2)`, and optional hyperprior as well.  Additionally it can be fixed to a constant by setting `nug.p = 0` and specifying the staring value as `gd[1]`.

* The mean is not fixed at zero but is instead a linear model, we need hyperparameters for it as well.  The key parameter there is the slope $\beta$.  The functional form for it's prior is specified as `bprior`, which takes values `bflat`, `b0`, `bmzt`, or `bmznot`.  A sensible choice is as a Gaussian `b0`, as shown in Equation 1 of the vignette

$$
\begin{align}
Z | \beta, \sigma^2, K &\sim N_n(\mathbf{F} \beta, \sigma^2 \mathbf{K}) \\
\beta | \sigma^2, \tau^2, \mathbf{W}, \beta_0 &\sim N_m(\beta_0, \sigma^2 \tau^2 \mathbf{W} )  \\
\beta_0 &\sim N_m(\mathbf{\mu}, \mathbf{B})
\end{align}
$$

## Example comparison

Say, given some observed data pairs $x,y$ and predictions desired on a grid $X$,

```{r}
X = seq(-5, 5, length=50)
x = c(-4, -3, -1,  0,  2)
y = c(-2,  0,  1,  2, -1)
```

Manually I could just do:

```{r}
d = .5; epsilon = .1; sigma = 1; #fixed hyperparamaters
SE <- function(Xi,Xj, d) sigma * exp(-0.5 * (Xi - Xj) ^ 2 / d ^ 2)
cov <- function(X, Y) outer(X, Y, SE, d)
cov_xx_inv <- solve(cov(x, x) + epsilon * diag(1, length(x)))
Ef <- cov(X, x) %*% cov_xx_inv %*% y
Cf <- cov(X, X) - cov(X, x)  %*% cov_xx_inv %*% cov(x, X)
```

And `Ef` and `Cf` are the estimated (krieging) mean and covariance.

## Priors

Note: it appears that the mixed gammas for `d` and `nug` actually treat the second argument, `g1` and `g2`, as the rate instead of the scale, so I invert the desired scale:

```{r}
s2.p <- c(50,50)
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)
```

Define as curves and plot the priors:

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

xx <- seq(.0001, 5, length.out=100)
priors <- data.frame(x = xx, s2 = s2_prior(xx), tau2 = tau2_prior(xx), beta0 = beta0_prior(xx, 1), nug = nug_prior(xx), d = d_prior(xx))
priors <- melt(priors, id="x")
ggplot(priors) + geom_line(aes(x, value)) + facet_wrap(~variable, scale="free")
```

```{r}
gp <- bgp(X=x, XX=X, Z=y, verb=0,
          meanfn="constant", bprior="b0", BTE=c(2000,6000,2), m0r1=FALSE,
          corr="exp", trace=TRUE, beta = 0,
          s2.p = s2.p, d.p = d.p, nug.p = nug.p,
          s2.lam = "fixed", d.lam = "fixed", nug.lam = "fixed",
          tau2.lam = "fixed", tau2.p = tau2.p)
```

Extract the posterior Gaussian process mean and the $\pm 2$ standard deviations over the predicted grid from the fit:

```{r}
V <- gp$ZZ.ks2
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))
```

and likewise for our manual data:

```{r}
manual_dat <- data.frame(x = X,
                         y = Ef,
                         ymin = (Ef - 1.96 * sqrt(diag(Cf))),
                         ymax = (Ef + 1.96 * sqrt(diag(Cf))))
```

Compare the GP posteriors:

```{r gp-plot}
ggplot(tgp_dat) +
    geom_ribbon(aes(x, y, ymin = ymin, ymax = ymax), fill="red", alpha = .1) + # Var
    geom_line(aes(x, y), col="red") + # mean
    geom_point(data = data.frame(x = x, y = y), aes(x, y)) + # raw data
    geom_ribbon(dat = manual_dat, aes(x, y, ymin = ymin, ymax = ymax), fill = "blue", alpha = .1) + # Var
    geom_line(dat = manual_dat, aes(x, y), col = "blue")  + #MEAN
    theme_bw() + theme(plot.background = element_rect(fill = "transparent",colour = NA))
```

## Posteriors of hyperparameters

```{r}
require(reshape2)
hyperparameters <- c("index", "s2", "tau2", "beta0", "nug", "d", "ldetK")
posteriors <- melt(gp$trace$XX[[1]][,hyperparameters], id="index")
ggplot(posteriors) + geom_histogram(aes(value)) + facet_wrap(~variable, scales="free")

```


text:
# Learning ADMB

Plotting and knitr options, (can generally be ignored)

### Model and parameters

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

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

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

### Sample Data

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

![plot of chunk obs](http://farm3.staticflickr.com/2807/8956942302_0d7d47ea49_o.png)

## Maximum Likelihood "by hand"

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

  f
  }
```

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

```r
init <- c(p, sigma_g)
mloglik(init) #true minus loglik
```

```
[1] -35.72
```

```r
o <- optim(init, mloglik, method="L", lower=1e-5, upper=1e5)
o$value
```

```
[1] -247.6
```

```r
o$par
```

```
[1] 0.9967297 9.9813554 5.1742699 0.0008183
```

While starting from arbitrary values we still find the optim.

```r
init <- c(1,1,1,1)
init <- c(p, sigma_g)
mloglik(init) #true minus loglik
```

```
[1] -35.72
```

```r
o <- optim(init, mloglik, method="L", lower=1e-5, upper=1e5)
o$value
```

```
[1] -247.6
```

```r
o$par
```

```
[1] 0.9967297 9.9813554 5.1742699 0.0008183
```

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

```r
# install_github("R2admb", "bbolker", subdir="R2admb") # dev version
library(R2admb)
```

## ADMB definition

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

The equivalent model

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

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

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

```r
setup_admb("/var/admb")
```

```
[1] "/var/admb"
```

```r

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

```
Model file: model_gen
Negative log-likelihood: 146.5
Coefficients:
     r      k      c      s
0.9992 1.0023 1.0004 9.4369
```

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

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

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

```
Warning: running command './model_gen > model_gen.out' had status 1
```

```r
m1
```

```
Model file: model_gen
Negative log-likelihood: -423.8
Coefficients:
        r         k         c         s
2.025e-10 9.498e+00 1.051e+01 1.000e-05
```

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

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

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


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

We use a Beverton-Holt model to drive the underlying dynamics, with parameters $A =$ `1.5` and $B =$ `0.05`.

```r
sigma_g <- 0.02
z_g <- function(sigma_g) rlnorm(1, 0, sigma_g) #1+(2*runif(1, 0,  1)-1)*sigma_g #
x_grid <- seq(0, 1.5 * K, length=101)
h_grid <- x_grid
profit = function(x,h) pmin(x, h)
delta <- 0.01
OptTime = 20
xT = 0
reward = profit(x_grid[length(x_grid)], x_grid[length(x_grid)]) + 1 / (1 - delta) ^ OptTime
## 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
x_0_observed <- x_grid[2]
```

We consider stochastic growth driven by a lognormal noise process, $X_{t+1} = z_g f(X_t)$, where $f$ is the stock recruitment curve and $z_g$ a lognormal shock with $\sigma_g$ = `0.02`.

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

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

Simulate data

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

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

```r
s2.p <- c(50,50)
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)
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)
```

Estimates a Ricker curve with parameters $r =$ `0.4182` and $K =$ `7.7853`

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

```
Warning: for memory/storage reasons, trace not recommended when
3*(10+d)*(BTE[2]-BTE[1])*R*(nn+1)/BTE[3]=23562000 > 1e+7.  Try reducing
nrow(XX)
```

We fit a Gaussian process with

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

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

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

```r
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 = dnorm, 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")
```

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

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

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

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

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

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

![plot of chunk gp-posteriors](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-gp-posteriors1.png)

```r
ggplot(prior_curves) +
  geom_line(aes(x=value, y=density), col="red", lwd=2) +
  facet_wrap(~ variable, scale="free")
```

![plot of chunk gp-posteriors](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-gp-posteriors2.png)

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

## The transition matrix of the inferred process

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

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

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

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

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

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

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

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

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

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

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

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

```r
policies <- melt(data.frame(stock=x_grid,
                            GP = x_grid[opt_gp$D[,1]],
                            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
```

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

```r
z_g = function() rlnorm(1, 0, sigma_g)
z_m = function() 1+(2*runif(1, 0,  1)-1) * sigma_m
```

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

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

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

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

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

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

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

```r
profits <- dt[, sum(profit), by = c("reps", "method")]
means <- profits[, mean(V1), by = method]
sds <- profits[, sd(V1), by = method]
cbind(means, sd = sds$V1)
```

```
       method    V1     sd
1:         GP 23.83 0.3527
2: Parametric 23.65 0.4774
3:       True 23.68 0.4513
```

Myers RA, Barrowman NJ, Hutchings JA and Rosenberg AA (1995).
“Population Dynamics of Exploited Fish Stocks at Low Population Levels.”
Science, 269.
ISSN 0036-8075, http://dx.doi.org/10.1126/science.269.5227.1106.