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

We use the model of Myers _et. al._ (1995).

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

With parameters `1, 2, 4.5`.

```r
x_0_observed <- allee + x_grid[5]
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-01d0f47bcc-sim-obs.png)

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

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

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

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

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

```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-01d0f47bcc-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") +
  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-01d0f47bcc-gp-posteriors1.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)
```

```
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-01d0f47bcc-gp-posteriors2.png)

## The transition matrix of the inferred process

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

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

![plot of chunk gp-F-sim](http://carlboettiger.info/assets/figures/2012-12-19-01d0f47bcc-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-01d0f47bcc-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-01d0f47bcc-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-01d0f47bcc-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-01d0f47bcc-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-01d0f47bcc-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 2.992 0.4489
2: Parametric 0.000 0.0000
3:       True 3.114 0.4534
```

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:
## Introductory examples in approximate dynamic programming

_Based on Powell 2006, page 97.  I try to conform to that notation throughout_

## Setup

First we define the Beverton-Holt stock-recruitment relationship as a function of stock size `x`, harvest `h` and parameters `p`

```r
f <- function(x, h, p){
    A <- p[1]
    B <- p[2]
    s <- pmax(x-h, 0)
    A * s/(1 + B * s)
}
p <- pars <- c(1.5, 0.5)
K <- (p[1] - 1)/p[2]
sigma_g <- 0.2
```

We begin with a simulation method $X_{t+1} = f(X_t, Z_t)$.  For illustration, let us consider $f(X_t, Z_t) = Z_t \frac{a X_t}{b + X_t}$ with a = 1.5 and b = 0.5.  We define a statespace $S$

```r
S <- seq(0, 1, length=11)
```

as a uniform grid of 11 points from 0 to 1.
We also need a value function on the state space, $C_t(S_t)$.
For simplicity, we set the price of harvest at unity and
the cost of harvesting at zero, so that $C_t(S_t, x_t) = \min(x_t, S_t)$.
($C_t$ is sometimes denoted $\mathbb{\Pi}$).
We also need an action space $\chi_t$ of possible harvest values.
Again for simplicity we assume that harvest can be set to any possible state size, $\chi_t \equiv S_t$,

```r
chi <- S
```

```r
T <- 10
N <- 10
```

The approximate dynamic programming algorithm will perform a finite number $N$ = 10 iterations over a window of time $T$ =10 in our example.  The algorithm can then be described as follows:

## Algorithm (1)

- **Step 0**

  - Initialize some value $\tilde V_t^0(S_t)$ for all states $S_t$,
  where the superscripts denote iterations in the forward approximation.
  As we know absolutely nothing yet to base our initial guess on,
  we just arbitrarily set this to zero.

```r
V <- numeric(length(S))
```

  - Choose some initial state $S_0^1$
  We start at some initial state for $n = 1$ (superscript) and $t = 0$ (subscript).
  The choice of initial condition may come from the problem itself,
  otherwise we choose something arbitrarily.

```r
S_0 <- 0.5
```

  - Set $n = 1$

- **Step 1**: Choose a sample path, $\omega^n$ (a vector of random draws)

```r
sigma <- 0.2
omega_n <- rlnorm(T, 0, sigma)
```

- ** Step 2**: For $t = 0, 1, 2, \ldots, T$, do:

  - Solve:

$$V_t(S_t) = \max_{x_t \in \chi_t} \left(C(S_t, x_t) + \gamma \sum_{s^{\prime} \in \mathcal{S}} \mathbb{P}(s^{\prime} | S_t^n, x_t) V_{t+1}^{n-1} s^{\prime} \right)$$

That is, choose action $x_t$ that maximizes the value of the next step.

Let's start with $t=0$, $n=1$ and fix an $x_0$ from the set of $\chi$
(allowing the action space to be the same in each period, we can omit
the subscript on $\chi$) to get started.  We first compute $C(S_0, x_0)$.

$S_0 = S_0^1$ which we fixed in step **0b** arbitrarily at 0.5.

The profits/costs $C(S_t, x_t)$ are the value derived by action (harvest) $x_t$
at state (stock) $S_t$.  Assuming a fixed price and no costs to harvesting,
this is just whichever number is smaller (since we cannot harvest more than
the available stock,

```r
C <- function(S, X) pmin(S, X)
```

(where we have used R's vectorized form of the min function).

This forward dynamic programming will still rely on the one-step transition
matrix, $\mathbb{P}$.

Let's get the trasition matrices for this problem, assuming log-normal noise,

```r
sdp_matrix <- determine_SDP_matrix(f, p, x_grid=S, h_grid=chi, sigma_g)
```

Which is a list of matrices, one for each harvest (action) $x_t$.

Then we want to consider a fixed $S_t^n$ and fixed $x_t$, and take the sum of
$\mathbb{P}(s^{\prime} | S_t^n, x_t)$ over the $s^{\prime}$, which means we want
the $x_t$ element from the list, and then we need sum over the distribution of future
states given the current state $S_t^n$, e.g. a row of the matrix, e.g. `sdp_matrix[[x]][s,]`,
which we (vector) multiply by $V_{t+1}^{n-1}(s^{\prime})$.

This value $V$ is of course unknown, other than our initial random guess $V_{t}^0$.
As we step through the iterations $V_t^1$, $V_t^2$, $V_t^3$, etc., this should convgerge to
something meaningful.

Note that the index along $S$ corresponding to $S_t^n$ is given by

```r
s <- which.min(abs(S-S_0))
```

So our maximization across $x$ just involves:

```r
values <-
  sapply(1:length(chi), function(x)
    C(S[s], chi[x]) + sdp_matrix[[x]][s,] %*% V
)

max_x <- which.max(values)
v_hat <- max(values)
```

Trivially, this is just the harvest level that maximizes $C$ so far (which
is just harvesting the $S_0$, since $\bar V^0_t$ begins at zero:

```r
chi[max_x]
```

```
[1] 0.5
```

- step **2b** We can now update our $\bar V^0_t$ to get $\bar V^1_t$, using the rule:

$$V_t^n(S_t) = \begin{cases}
\hat v_t^n & S_t = S_t^n \\
\bar V_t^{n-1}(S_t) & \textrm{otherwise}
\end{cases}$$

e.g. use our maximum value for the case of the state we just considered $S_t = S_t^n$, otherwise leave $V_t$ unchanged.  Our new $V$ is thus:

```r
V[s] = v_hat
```

- step **2c**  Compute $S^n_{t+1} = S^M(S_t^n, x^n_t, W_{t+1}(\omega^n))$

We compute the next state using our `max_x` for $x^n_t$, our random samples
and the transition function...

```r
S_1 <- omega_n[1] * f(S_0, chi[x_nt], p)
```

```
Error: object 'x_nt' not found
```

- **Step 3** Let $n = n+1$. if $n < N$, go to step 1

assembling the solution...

## As a recursive algorithm

Things we want to keep track of:

* $V^{n+1}_t$

We have time, iterations, and state/action space to keep track over.  eek.

```r
N <- 5000 # iterations
M <- 20 # gridsize
Tmax <- 5 # Time horizon

gamma <- 0.95 # Discount
# f
# p

sigma_g <- 0.5 # larger variation in random draws helps
chi <- seq(0, 1, length.out = M)
S <- seq(0, 1, length.out = M)

sdp_matrix <- determine_SDP_matrix(f, p, x_grid=S, h_grid=chi, sigma_g)

V <- matrix(1, M, Tmax)  # A* strategy
# Fails to explore at matrix(0, M, Tmax)
# consider: # V <- matrix(rep(chi, Tmax), nrow=M) #
# V[,1] <- chi   # fails to explore if it doesn't have at least some non-zero values

C <- function(S, X) pmin(S, X)
S_0 <- 0.5
alpha <- 1 # learning rate

for(n in 1:N){

  omega_n <- rlnorm(Tmax, 0, sigma_g)
  S_current <- S_0 #runif(1,0,1) # explores faster when this is random

  for(t in 1:Tmax){
    # index of the state we're considering
    s <- which.min(abs(S-S_current))

    # Find the action maximizing the value
    values <- sapply(1:length(chi), function(x)
      C(S[s], chi[x]) + gamma * sdp_matrix[[x]][s,] %*% V[,t])
    hat <-  c(x_nt = which.max(values), v_nt = max(values))

    # Update value V as mixture of new value and previous value
    V[hat["x_nt"], t] <- (1 - alpha) * V[hat["x_nt"],t] + alpha * hat["v_nt"]

    # Advance the state in time along random path
    S_current <- omega_n[t] * f(S_current, chi[hat["x_nt"]], p)

  }
}
```

for comparison: the SDP solution

```r
opt <- find_dp_optim(sdp_matrix, S, chi, 70, 0.5, C, 1-gamma, reward=0)
opt$V
```

```
 [1] 0.000 7.173 7.375 7.496 7.584 7.653 7.710 7.760 7.810 7.860 7.910
[12] 7.960 8.010 8.060 8.110 8.160 8.210 8.260 8.310 8.360
```

### Problems arising from the discretization

Note that after the first iteration, $n=1$, the value matrix $V$ is no longer all zeros.  There is a single state, $S = S_0 =$ 0.5, at which we have value.  That value is lost if we set harvest $x$ too high, since we know we will not then end up in that state -- from whence comes the incentive to consider future value.  Unfortunately, the value exists only if we hit that state exactly -- all other states are assumed to have zero value still.

### Additional problems

We no longer have the loop-over-all-states problem, but we face several new or remaining issues:

1. We still require the use of the one-step transition matrix, with the equally troublesome sum over all states $\sum_{s^{\prime}\in S} \mathbb{P}(s^{\prime} | S_t^n, x_t)$.  We will fix this by approximating the transitions in step 2b using random draws as well.

2. We only update the values of states we visit.  We still need a way to estimate the value of states we have not visited.

3. Worse, we might not visit states that seem bad relative to states we have visited. This is particularly atrocious in this example.  Since we initialize the value of all states at 0, the algorithm prefers to harvest all stock from the current state rather than risk a transition into a state starting at 0. There is no convergence guarentee that we will ever escape this cycle of avoiding states we have not seen. We can alter the initial guess of the value of course, and we could alter the starting condition to better explore.

### Using Non-stochastic transition information only, step **2b** can be written as:

Taking $x_0$ as the smallest harvest, $\min(\chi)$ = 0 and
evaluating $C(S_0,X_0) = \min(S_0, X_0)$ gives us 0,
rather trivially.
The next terms depend on the value $\tilde V^0_1(s^{\prime})$ for
all $s^{\prime} \in S$, which we have no idea about.  Fortunately we have
assumed a value for each of these in step 0a.

We must also come up with some values for the probability
$\mathbb{P}(s^{\prime} | S_1^0, x_1)$ for each state, given our current
state $S_1^0$ and considered action $x_1$.  This is more straight forward,
since it is determined by our one-step transition function (without
simulation - recall that the single step transition is given exactly).

To do so, we evaluate the argument for each value in our action
space, $x_t \in \chi_t$,

```r
s <- S_0
C <- function(S, X) pmin(S, X)
arg <- sapply(chi, function(x) C(s, x) + f(S, x, p) %*% V)
x_nt = which.max(arg)
v_nt = max(arg)
V[x_nt] = v_nt
```

repo_name: cboettig/nonparametric-bayes

repo_language: TeX

mime_type: text/plain


text:
## Optimization over hyperparameters:

```{r results="hide"}
opts_knit$set(upload.fun = socialR::notebook.url)
```

Static definition of gp function

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

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

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

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

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

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

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

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

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

Example plot

```{r}
plot.gpfit(gp_fit(obs, X, pars=c("sigma_n"=.3, "l"=.6)))
```

estimates non-boundary solution of `l`:

```{r}
minusloglik <- function(par){
  gp <- gp_fit(obs, X, c(sigma_n=.3, l=par))
  -gp$llik
}
minusloglik(1)
o <- optimize(minusloglik, c(0,40))
o

```

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

```{r}
minusloglik <- function(par){
  gp <- gp_fit(obs, X, c(sigma_n=par, l=.6))
  -gp$llik
}
minusloglik(1)
o <- optimize(minusloglik, c(0,40))
o

```

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

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

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

Weirdly inconsistent?

```{r}
gp_k <- gausspr(obs$x, obs$y, kernel="rbfdot")
1/sqrt(2*as.numeric(gp_k@kernelf@kpar))
gp_k <- gausspr(obs$x, obs$y, kernel="rbfdot")
1/sqrt(2*as.numeric(gp_k@kernelf@kpar))
```

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

```{r results="hide"}
knit("beverton_holt_data.Rmd")
```

```{r}
minusloglik <- function(par){
  if(any(par < 0))
    Inf
  else {
  gp <- gp_fit(obs, X, par)
  -gp$llik
  }
}
par <- c(sigma_n=1, l=1)
minusloglik(par)
#o <- optim(par, minusloglik)
o <- optim(par, minusloglik, method="L", lower=c(0,0), upper=c(20,20))
o
hyperpars <- o$par
```

**Yikes.**

let's try 1-D optimization:

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

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

Plot the optimal hyperparameter solution:

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


text:
```r
## Set the paths for cache and figure
library(methods)
library(knitr)
basename <- gsub(".Rmd", "", knitr:::knit_concord$get("infile"))
opts_chunk$set(fig.path = paste("figure/", basename, "-", sep = ""), cache.path = paste("cache/",
    basename, "/", sep = ""))
opts_chunk$set(cache = 1)
opts_chunk$set(tidy = FALSE, warning = FALSE, message = FALSE, comment = NA,
    verbose = TRUE)
```

```r
library(nonparametricbayes)

models <- c("Myers","Allen")

parameters <- list(Myers = list(
  c(r=1.5 + rnorm(1, 0, .1), theta=2.5 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .2)),
  c(r=1.5 + rnorm(1, 0, .1), theta=2.5 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .2)),
  c(r=1.5 + rnorm(1, 0, .1), theta=2.5 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .2))),
                   Allen = list(
  c(r=2 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .1), C=5 + rnorm(1, 0, .2)),
  c(r=2 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .1), C=5 + rnorm(1, 0, .2)),
  c(r=2 + rnorm(1, 0, .1), K=8 + rnorm(1, 0, .1), C=5 + rnorm(1, 0, .2)))
                   )
nuisance_pars <- c("sigma_g")
nuisance_values <- list(sigma_g = c(0.01, 0.05, 0.1))
replicates <- c(1111, 2222, 3333, 4444, 5555, 6666, 7777, 8888) # seeds

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

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

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

  sigma_g <- nuisance[["sigma_g"]]
  z_g <- function() rlnorm(1, 0, sigma_g)
  p <- parameters

    f = f
    p = p

    z_g = z_g
    x_grid = seq(0, 15, length=50)
    h_grid = x_grid
    profit = function(x,h) pmin(x, h)
    delta = 0.01
    OptTime = 50  # stationarity with unstable models is tricky thing
    reward = 0
    xT = 0
    Xo = 5.5 # observations start from
    x0 = 8 # simulation under policy starts from
    Tobs = 40
    MaxT = 1000 # timeout for value iteration convergence

  # replicate over random seed
  yields <- sapply(seed,
                   function(seed_i){
                     set.seed(seed_i)

## Simulate data
                     x <- numeric(Tobs)
                     x[1] <- Xo
                     nz <- 1
                     for(t in 1:(Tobs-1))
                       x[t+1] = z_g() * f(x[t], h=0, p=p)

                     X = c(rep(0,nz), pmax(rep(0,Tobs-1), x[1:(Tobs-1)]))
                     Y = c(rep(0,nz), x[2:Tobs])

## GP Stuff

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

                     ## @knitr gp
                     gp <- gp_mcmc(X, y=Y, n=1e5, s2.p = s2.p, d.p = d.p)
                     gp_dat <- gp_predict(gp, x_grid, burnin=1e4, thin=300)

                     matrices_gp <- gp_transition_matrix(gp_dat$Ef_posterior, gp_dat$Vf_posterior, x_grid, h_grid)
                     opt_gp <- value_iteration(matrices_gp, x_grid, h_grid, MaxT, xT, profit, delta, reward)

## @knitr jags-setup
y <- x
N <- length(x);
jags.data <- list("N","y")
n.chains <- 6
n.iter <- 1e6
n.burnin <- floor(10000)
n.thin <- max(1, floor(n.chains * (n.iter - n.burnin)/1000))
n.update <- 10

## @knitr common-priors
stdQ_prior_p <- c(1e-6, 100)
stdR_prior_p <- c(1e-6, .1)
stdQ_prior  <- function(x) dunif(x, stdQ_prior_p[1], stdQ_prior_p[2])
stdR_prior  <- function(x) dunif(x, stdR_prior_p[1], stdR_prior_p[2])

## @knitr allen-model
K_prior_p <- c(0.01, 20.0)
r0_prior_p <- c(0.01, 6.0)
theta_prior_p <- c(0.01, 20.0)

bugs.model <-
paste(sprintf(
"model{
  K     ~ dunif(%s, %s)
  r0    ~ dunif(%s, %s)
  theta ~ dunif(%s, %s)
  stdQ ~ dunif(%s, %s)",
  K_prior_p[1], K_prior_p[2],
  r0_prior_p[1], r0_prior_p[2],
  theta_prior_p[1], theta_prior_p[2],
  stdQ_prior_p[1], stdQ_prior_p[2]),

  "
  iQ <- 1 / (stdQ * stdQ);
  y[1] ~ dunif(0, 10)
  for(t in 1:(N-1)){
    mu[t] <- log(y[t]) + r0 * (1 - y[t]/K)* (y[t] - theta) / K
    y[t+1] ~ dlnorm(mu[t], iQ)
  }
}")
writeLines(bugs.model, "allen_process.bugs")

## @knitr allen-priors
K_prior     <- function(x) dunif(x, K_prior_p[1], K_prior_p[2])
r0_prior <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2])
theta_prior <- function(x) dunif(x, theta_prior_p[1], theta_prior_p[2])
par_priors  <- list(K = K_prior, deviance = function(x) 0 * x,
                    r0 = r0_prior, theta = theta_prior,
                    stdQ = stdQ_prior)

## @knitr allen-mcmc
jags.params=c("K","r0","theta","stdQ") # be sensible about the order here
jags.inits <- function(){
  list("K"= 10 * rlnorm(1,0, 0.1),
       "r0"= 1 * rlnorm(1,0, 0.1) ,
       "theta"=   5 * rlnorm(1,0, 0.1) ,
       "stdQ"= abs( 0.1 * rlnorm(1,0, 0.1)),
       .RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}

set.seed(1234)
# parallel refuses to take variables as arguments (e.g. n.iter = 1e5 works, but n.iter = n doesn't)
allen_jags <- do.call(jags, list(data=jags.data, inits=jags.inits,
                                      jags.params, n.chains=n.chains,
                                      n.iter=n.iter, n.thin=n.thin,
                                      n.burnin=n.burnin,
                                      model.file="allen_process.bugs"))

# Run again iteratively if we haven't met the Gelman-Rubin convergence criterion
recompile(allen_jags) # required for parallel
allen_jags <- do.call(autojags,
											list(object=allen_jags, n.update=n.update,
                           n.iter=n.iter, n.thin = n.thin))

## @knitr allen-traces
tmp <- lapply(as.mcmc(allen_jags), as.matrix) # strip classes the hard way...
allen_posteriors <- melt(tmp, id = colnames(tmp[[1]]))
## @knitr allen-output
A <- allen_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
pardist <- acast(A, index ~ variable)

## @knitr ricker-model
K_prior_p <- c(0.01, 40.0)
r0_prior_p <- c(0.01, 20.0)
bugs.model <-
paste(sprintf(
"model{
  K    ~ dunif(%s, %s)
  r0    ~ dunif(%s, %s)
  stdQ ~ dunif(%s, %s)",
  K_prior_p[1], K_prior_p[2],
  r0_prior_p[1], r0_prior_p[2],
  stdQ_prior_p[1], stdQ_prior_p[2]),
  "
  iQ <- 1 / (stdQ * stdQ);
  y[1] ~ dunif(0, 10)
  for(t in 1:(N-1)){
    mu[t] <- log(y[t]) + r0 * (1 - y[t]/K)
    y[t+1] ~ dlnorm(mu[t], iQ)
  }
}")
writeLines(bugs.model, "ricker_process.bugs")

## @knitr ricker-priors
K_prior     <- function(x) dunif(x, K_prior_p[1], K_prior_p[2])
r0_prior <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2])
par_priors <- list(K = K_prior, deviance = function(x) 0 * x,
                   r0 = r0_prior, stdQ = stdQ_prior)

## @knitr ricker-mcmc
jags.params=c("K","r0", "stdQ")
jags.inits <- function(){
  list("K"= 10 * rlnorm(1,0,.5),
       "r0"= rlnorm(1,0,.5),
       "stdQ"=sqrt(0.05) * rlnorm(1,0,.5),
       .RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
set.seed(12345)
ricker_jags <- do.call(jags,
                       list(data=jags.data, inits=jags.inits,
                            jags.params, n.chains=n.chains,
                            n.iter=n.iter, n.thin=n.thin, n.burnin=n.burnin,
                            model.file="ricker_process.bugs"))
recompile(ricker_jags)
ricker_jags <- do.call(autojags,
                       list(object=ricker_jags, n.update=n.update,
														n.iter=n.iter, n.thin = n.thin,
														progress.bar="none"))

## @knitr ricker-traces
tmp <- lapply(as.mcmc(ricker_jags), as.matrix) # strip classes the hard way...
ricker_posteriors <- melt(tmp, id = colnames(tmp[[1]]))
names(ricker_posteriors) = c("index", "variable", "value", "chain")

## @knitr ricker-output
A <- ricker_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
ricker_pardist <- acast(A, index ~ variable)

## @knitr myers-model
r0_prior_p <- c(.0001, 10.0)
theta_prior_p <- c(.0001, 10.0)
K_prior_p <- c(.0001, 40.0)
bugs.model <-
paste(sprintf(
"model{
  r0    ~ dunif(%s, %s)
  theta    ~ dunif(%s, %s)
  K    ~ dunif(%s, %s)
  stdQ ~ dunif(%s, %s)",
  r0_prior_p[1], r0_prior_p[2],
  theta_prior_p[1], theta_prior_p[2],
  K_prior_p[1], K_prior_p[2],
  stdQ_prior_p[1], stdQ_prior_p[2]),

  "
  iQ <- 1 / (stdQ * stdQ);

  y[1] ~ dunif(0, 10)
  for(t in 1:(N-1)){
    mu[t] <- log(r0)  + theta * log(y[t]) - log(1 + pow(abs(y[t]), theta) / K)
    y[t+1] ~ dlnorm(mu[t], iQ)
  }
}")
writeLines(bugs.model, "myers_process.bugs")

## @knitr myers-priors
K_prior     <- function(x) dunif(x, K_prior_p[1], K_prior_p[2])
r_prior     <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2])
theta_prior <- function(x) dunif(x, theta_prior_p[1], theta_prior_p[2])
par_priors <- list( deviance = function(x) 0 * x, K = K_prior,
                    r0 = r_prior, theta = theta_prior,
                    stdQ = stdQ_prior)

## @knitr myers-mcmc
jags.params=c("r0", "theta", "K", "stdQ")
jags.inits <- function(){
  list("r0"= 1 * rlnorm(1,0,.1),
       "K"=    10 * rlnorm(1,0,.1),
       "theta" = 1 * rlnorm(1,0,.1),
       "stdQ"= sqrt(0.2) * rlnorm(1,0,.1),
       .RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
set.seed(12345)
myers_jags <- do.call(jags,
                      list(data=jags.data, inits=jags.inits,
													 jags.params, n.chains=n.chains,
													 n.iter=n.iter, n.thin=n.thin,
                           n.burnin=n.burnin,
                           model.file="myers_process.bugs"))
recompile(myers_jags)
myers_jags <- do.call(autojags,
                      list(myers_jags, n.update=n.update,
                           n.iter=n.iter, n.thin = n.thin,
                           progress.bar="none"))

## @knitr myers-traces
tmp <- lapply(as.mcmc(myers_jags), as.matrix) # strip classes
myers_posteriors <- melt(tmp, id = colnames(tmp[[1]]))
names(myers_posteriors) = c("index", "variable", "value", "chain")
## @knitr myers-output
A <- myers_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
myers_pardist <- acast(A, index ~ variable)

allen_deviance  <- - posterior.mode(pardist[,'deviance'])
ricker_deviance <- - posterior.mode(ricker_pardist[,'deviance'])
myers_deviance  <- - posterior.mode(myers_pardist[,'deviance'])
true_deviance   <- 2*estf(c(p, sigma_g))
mle_deviance    <- 2*estf(c(est$p, est$sigma_g))
aictable <- data.frame(Allen = allen_deviance + 2*(1+length(bayes_pars)),  # +1 for noise parameter
                       Ricker = ricker_deviance + 2*(1+length(ricker_bayes_pars)),
                       Myers = myers_deviance + 2*(1+length(myers_bayes_pars)),
                       row.names = c("AIC"))
bictable <- data.frame(Allen = allen_deviance + log(length(x))*(1+length(bayes_pars)),
                       Ricker = ricker_deviance + log(length(x))*(1+length(ricker_bayes_pars)),
                       Myers = myers_deviance + log(length(x))*(1+length(myers_bayes_pars)),
                       row.names = c("BIC"))

df <- rbind(dictable, aictable, bictable)

                   })

  dat <- melt(yields, id=names(yields[[1]]))
}
```

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

```
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 328

Initializing model

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

Initializing model
```

```
Error: replacement has 0 rows, data has 5070
```

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

```
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 328

Initializing model

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

Initializing model
```

```
Error: replacement has 0 rows, data has 5070
```

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

```
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 328

Initializing model

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

Initializing model
```

```
Error: replacement has 0 rows, data has 5070
```

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

```
Error: system is computationally singular: reciprocal condition number =
1.82252e-16
```

```
Timing stopped at: 0.144 0.004 0.148
```

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

```
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 328

Initializing model

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

Initializing model
```

```
Error: replacement has 0 rows, data has 5070
```

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

```
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 328

Initializing model

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

Initializing model
```

```
Error: replacement has 0 rows, data has 5070
```

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

```
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 328

Initializing model

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

Initializing model
```

```
Error: replacement has 0 rows, data has 5070
```

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

```
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 328

Initializing model

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

Initializing model
```

```
Error: replacement has 0 rows, data has 5070
```

```r
## Assemble into data.frame
allen_dat <- rbind(allen1.01, allen1.05,
             allen2.01, allen2.05)
```

```
Error: object 'allen1.01' not found
```

```r
myers_dat <- rbind(Myers1.01, Myers1.05,
             Myers2.01, Myers2.05)
```

```
Error: object 'Myers1.01' not found
```


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

```{r libraries, include=FALSE}
```

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

```{r BevHolt}
```

We use a Beverton-Holt model to drive the underlying dynamics, with parameters $A =$ `r p[1]` and $B =$ `r p[2]`.

```{r sdp-pars-bh}
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$ = `r sigma_g`.

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

Simulate data

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

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

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

Estimates a Ricker curve with parameters $r =$ `r p_alt[1]` and $K =$ `r p_alt[2]`

```{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 est-F-sim, dependson="persistence-test"}
```

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

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

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

```{r policy_plot}
```

```{r othernoise}
```

```{r stationary_policy_only}
```

```{r simulate}
```

```{r tidy}
```

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

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

```{r costs}
```

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


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

![plot of chunk unnamed-chunk-6](http://carlboettiger.info/assets/figures/2012-11-19-58073767b5-unnamed-chunk-6.png)


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

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

```
       x variable   value
1 -5.000       X1 -0.4055
2 -4.796       X1 -0.6389
3 -4.592       X1 -0.8300
4 -4.388       X1 -0.9550
5 -4.184       X1 -0.9996
6 -3.980       X1 -0.9581
```

Plot the result

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

![plot of chunk fig2a](http://farm8.staticflickr.com/7283/8738379037_e1314a3c34_o.jpg)

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

![plot of chunk fig2b](http://farm8.staticflickr.com/7288/8738379253_b298f9e11a_o.jpg)

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

![plot of chunk fig2_data](http://farm8.staticflickr.com/7281/8739499454_0b63134241_o.jpg)

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

![plot of chunk fig2c](http://farm8.staticflickr.com/7284/8739499586_dfe825b0e4_o.jpg)

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

```
       l  sigma.n
0.703891 0.002685
```

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

```
       [,1]
[1,] -21.87
```

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

```
       [,1]
[1,] -40.78
```

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

```
[1] 0.1509
```

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

![plot of chunk gp_mcmc](http://farm8.staticflickr.com/7285/8738380125_14bb7a3296_o.jpg)

```r

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

![plot of chunk gp_mcmc](http://farm8.staticflickr.com/7285/8739500886_70b018f8ab_o.jpg)

```r

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

![plot of chunk gp_posteriors](http://farm8.staticflickr.com/7281/8738381585_dc0bd82493_o.jpg)


text:
## Introductory examples in approximate dynamic programming

```{r  setup, echo=FALSE, cache=FALSE, results="hide", message=FALSE, warning=FALSE}
require(multipleuncertainty)
require(pdgControl)
require(reshape2)
require(ggplot2)
require(data.table)
````

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

opts_chunk$set(fig.path = paste("figure/",
                                gsub(".Rmd", "",
                                     knitr:::knit_concord$get('infile')),
                                "-", sep=""),
               cache.path = paste(gsub(".Rmd", "",
                                       knitr:::knit_concord$get('infile') ),
                                "/", sep="")
              )
require(ggplot2)
theme_notebook <-
  theme_grey() +
  theme(plot.background = element_rect(fill = "transparent", colour = NA),
        panel.grid.minor = element_line(colour = "transparent"))
theme_set(theme_notebook)
```

_Based on Powell 2006, page 97.  I try to conform to that notation throughout_

## Setup

First we define the Beverton-Holt stock-recruitment relationship as a function of stock size `x`, harvest `h` and parameters `p`

``` {r bevholt}
f <- function(x, h, p){
    A <- p[1]
    B <- p[2]
    s <- pmax(x-h, 0)
    A * s/(1 + B * s)
}
p <- pars <- c(1.5, 0.5)
K <- (p[1] - 1)/p[2]
sigma_g <- 0.2
```

We begin with a simulation method $X_{t+1} = f(X_t, Z_t)$.  For illustration, let us consider $f(X_t, Z_t) = Z_t \frac{a X_t}{b + X_t}$ with a = `r pars[1]` and b = `r pars[2]`.  We define a statespace $S$

```{r}
S <- seq(0, 1, length=11)
```

as a uniform grid of `r length(S)` points from `r min(S)` to `r max(S)`.
We also need a value function on the state space, $C_t(S_t)$.
For simplicity, we set the price of harvest at unity and
the cost of harvesting at zero, so that $C_t(S_t, x_t) = \min(x_t, S_t)$.
($C_t$ is sometimes denoted $\mathbb{\Pi}$).
We also need an action space $\chi_t$ of possible harvest values.
Again for simplicity we assume that harvest can be set to any possible state size, $\chi_t \equiv S_t$,

```{r action-space}
chi <- S
```

```{r problem-size}
T <- 10
N <- 10
```

The approximate dynamic programming algorithm will perform a finite number $N$ = `r N` iterations over a window of time $T$ =`r T` in our example.  The algorithm can then be described as follows:

## Algorithm (1)

- **Step 0**

  - Initialize some value $\tilde V_t^0(S_t)$ for all states $S_t$,
  where the superscripts denote iterations in the forward approximation.
  As we know absolutely nothing yet to base our initial guess on,
  we just arbitrarily set this to zero.

```{r}
V <- numeric(length(S))
```

  - Choose some initial state $S_0^1$
  We start at some initial state for $n = 1$ (superscript) and $t = 0$ (subscript).
  The choice of initial condition may come from the problem itself,
  otherwise we choose something arbitrarily.

```{r}
S_0 <- 0.5
```

  - Set $n = 1$

- **Step 1**: Choose a sample path, $\omega^n$ (a vector of random draws)

```{r}
sigma <- 0.2
omega_n <- rlnorm(T, 0, sigma)
```

- ** Step 2**: For $t = 0, 1, 2, \ldots, T$, do:

  - Solve:

$$V_t(S_t) = \max_{x_t \in \chi_t} \left(C(S_t, x_t) + \gamma \sum_{s^{\prime} \in \mathcal{S}} \mathbb{P}(s^{\prime} | S_t^n, x_t) V_{t+1}^{n-1} s^{\prime} \right)$$

That is, choose action $x_t$ that maximizes the value of the next step.

Let's start with $t=0$, $n=1$ and fix an $x_0$ from the set of $\chi$
(allowing the action space to be the same in each period, we can omit
the subscript on $\chi$) to get started.  We first compute $C(S_0, x_0)$.

$S_0 = S_0^1$ which we fixed in step **0b** arbitrarily at `r S_0`.

The profits/costs $C(S_t, x_t)$ are the value derived by action (harvest) $x_t$
at state (stock) $S_t$.  Assuming a fixed price and no costs to harvesting,
this is just whichever number is smaller (since we cannot harvest more than
the available stock,

```{r}
C <- function(S, X) pmin(S, X)
```

(where we have used R's vectorized form of the min function).

This forward dynamic programming will still rely on the one-step transition
matrix, $\mathbb{P}$.

Let's get the trasition matrices for this problem, assuming log-normal noise,

```{r transition-matrix}
sdp_matrix <- determine_SDP_matrix(f, p, x_grid=S, h_grid=chi, sigma_g)
```

Which is a list of matrices, one for each harvest (action) $x_t$.

Then we want to consider a fixed $S_t^n$ and fixed $x_t$, and take the sum of
$\mathbb{P}(s^{\prime} | S_t^n, x_t)$ over the $s^{\prime}$, which means we want
the $x_t$ element from the list, and then we need sum over the distribution of future
states given the current state $S_t^n$, e.g. a row of the matrix, e.g. `sdp_matrix[[x]][s,]`,
which we (vector) multiply by $V_{t+1}^{n-1}(s^{\prime})$.

This value $V$ is of course unknown, other than our initial random guess $V_{t}^0$.
As we step through the iterations $V_t^1$, $V_t^2$, $V_t^3$, etc., this should convgerge to
something meaningful.

Note that the index along $S$ corresponding to $S_t^n$ is given by

```{r}
s <- which.min(abs(S-S_0))
```

So our maximization across $x$ just involves:

```{r}
values <-
  sapply(1:length(chi), function(x)
    C(S[s], chi[x]) + sdp_matrix[[x]][s,] %*% V
)

max_x <- which.max(values)
v_hat <- max(values)
```

Trivially, this is just the harvest level that maximizes $C$ so far (which
is just harvesting the $S_0$, since $\bar V^0_t$ begins at zero:

```{r}
chi[max_x]
```

- step **2b** We can now update our $\bar V^0_t$ to get $\bar V^1_t$, using the rule:

$$V_t^n(S_t) = \begin{cases}
\hat v_t^n & S_t = S_t^n \\
\bar V_t^{n-1}(S_t) & \textrm{otherwise}
\end{cases}$$

e.g. use our maximum value for the case of the state we just considered $S_t = S_t^n$, otherwise leave $V_t$ unchanged.  Our new $V$ is thus:

```{r}
V[s] = v_hat
```

- step **2c**  Compute $S^n_{t+1} = S^M(S_t^n, x^n_t, W_{t+1}(\omega^n))$

We compute the next state using our `max_x` for $x^n_t$, our random samples
and the transition function...

```{r}
S_1 <- omega_n[1] * f(S_0, chi[x_nt], p)
```

- **Step 3** Let $n = n+1$. if $n < N$, go to step 1

assembling the solution...

## As a recursive algorithm

Things we want to keep track of:

* $V^{n+1}_t$

We have time, iterations, and state/action space to keep track over.  eek.

```{r}
N <- 5000 # iterations
M <- 20 # gridsize
Tmax <- 5 # Time horizon

gamma <- 0.95 # Discount
# f
# p

sigma_g <- 0.5 # larger variation in random draws helps
chi <- seq(0, 1, length.out = M)
S <- seq(0, 1, length.out = M)

sdp_matrix <- determine_SDP_matrix(f, p, x_grid=S, h_grid=chi, sigma_g)

V <- matrix(1, M, Tmax)  # A* strategy
# Fails to explore at matrix(0, M, Tmax)
# consider: # V <- matrix(rep(chi, Tmax), nrow=M) #
# V[,1] <- chi   # fails to explore if it doesn't have at least some non-zero values

C <- function(S, X) pmin(S, X)
S_0 <- 0.5
alpha <- 1 # learning rate

for(n in 1:N){

  omega_n <- rlnorm(Tmax, 0, sigma_g)
  S_current <- S_0 #runif(1,0,1) # explores faster when this is random

  for(t in 1:Tmax){
    # index of the state we're considering
    s <- which.min(abs(S-S_current))

    # Find the action maximizing the value
    values <- sapply(1:length(chi), function(x)
      C(S[s], chi[x]) + gamma * sdp_matrix[[x]][s,] %*% V[,t])
    hat <-  c(x_nt = which.max(values), v_nt = max(values))

    # Update value V as mixture of new value and previous value
    V[hat["x_nt"], t] <- (1 - alpha) * V[hat["x_nt"],t] + alpha * hat["v_nt"]

    # Advance the state in time along random path
    S_current <- omega_n[t] * f(S_current, chi[hat["x_nt"]], p)

  }
}
```

for comparison: the SDP solution

```{r}
opt <- find_dp_optim(sdp_matrix, S, chi, 70, 0.5, C, 1-gamma, reward=0)
opt$V
```

### Problems arising from the discretization

Note that after the first iteration, $n=1$, the value matrix $V$ is no longer all zeros.  There is a single state, $S = S_0 =$ `r S_0`, at which we have value.  That value is lost if we set harvest $x$ too high, since we know we will not then end up in that state -- from whence comes the incentive to consider future value.  Unfortunately, the value exists only if we hit that state exactly -- all other states are assumed to have zero value still.

### Additional problems

We no longer have the loop-over-all-states problem, but we face several new or remaining issues:

1. We still require the use of the one-step transition matrix, with the equally troublesome sum over all states $\sum_{s^{\prime}\in S} \mathbb{P}(s^{\prime} | S_t^n, x_t)$.  We will fix this by approximating the transitions in step 2b using random draws as well.

2. We only update the values of states we visit.  We still need a way to estimate the value of states we have not visited.

3. Worse, we might not visit states that seem bad relative to states we have visited. This is particularly atrocious in this example.  Since we initialize the value of all states at 0, the algorithm prefers to harvest all stock from the current state rather than risk a transition into a state starting at 0. There is no convergence guarentee that we will ever escape this cycle of avoiding states we have not seen. We can alter the initial guess of the value of course, and we could alter the starting condition to better explore.

### Using Non-stochastic transition information only, step **2b** can be written as:

Taking $x_0$ as the smallest harvest, $\min(\chi)$ = `r min(chi)` and
evaluating $C(S_0,X_0) = \min(S_0, X_0)$ gives us `r min(S_0, min(chi))`,
rather trivially.
The next terms depend on the value $\tilde V^0_1(s^{\prime})$ for
all $s^{\prime} \in S$, which we have no idea about.  Fortunately we have
assumed a value for each of these in step 0a.

We must also come up with some values for the probability
$\mathbb{P}(s^{\prime} | S_1^0, x_1)$ for each state, given our current
state $S_1^0$ and considered action $x_1$.  This is more straight forward,
since it is determined by our one-step transition function (without
simulation - recall that the single step transition is given exactly).

To do so, we evaluate the argument for each value in our action
space, $x_t \in \chi_t$,

```r
s <- S_0
C <- function(S, X) pmin(S, X)
arg <- sapply(chi, function(x) C(s, x) + f(S, x, p) %*% V)
x_nt = which.max(arg)
v_nt = max(arg)
V[x_nt] = v_nt
```


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

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

```{r graphing-options, include=FALSE}
opts_knit$set(upload.fun = socialR::notebook.url)
opts_chunk$set(dev.args=list(bg="transparent"), comment=NA, tidy=FALSE)
theme_set(theme_bw())
theme_update(panel.background = element_rect(fill = "transparent",colour = NA),
             plot.background = element_rect(fill = "transparent",colour = NA))
```

$$\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")
```

```{r phase-space, cache=TRUE, cache.path="phasespace/"}
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")
```

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 delay-equations}
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)
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")
```

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

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

`kernlab` method


text:
```{r set-options, echo = FALSE, cache = FALSE, external = TRUE, include = FALSE}
opts_chunk$set(external = TRUE, cache = FALSE, cache.path = "allee/", warning=FALSE)
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::flickr.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}
#inv gamma has mean b / (a - 1) (assuming a>1) and variance b ^ 2 / ((a - 2) * (a - 1) ^ 2) (assuming a>2)
s2.p <- c(5,5)
tau2.p <- c(5,1)
d.p = c(10, 1/0.1, 10, 1/0.1)
nug.p = c(10, 1/0.1, 10, 1/0.1) # gamma mean
s2_prior <- function(x) dinvgamma(x, s2.p[1], s2.p[2])
tau2_prior <- function(x) dinvgamma(x, tau2.p[1], tau2.p[2])
d_prior <- function(x) dgamma(x, d.p[1], scale = d.p[2]) + dgamma(x, d.p[3], scale = d.p[4])
nug_prior <- function(x) dgamma(x, nug.p[1], scale = nug.p[2]) + dgamma(x, nug.p[3], scale = nug.p[4])
beta0_prior <- function(x, tau) dnorm(x, 0, tau)
beta = c(0)
priors <- list(s2 = s2_prior, tau2 = tau2_prior, beta0 = dnorm, nug = nug_prior, d = d_prior, ldetK = function(x) 0)
```

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

```

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

```{r sdp-pars-explore}
sigma_g <- 0.05
sigma_m <- 0.0
z_g = function() rlnorm(1, 0, sigma_g)
z_m = function() 1+(2*runif(1, 0,  1)-1) * sigma_m
x_grid <- seq(0, 1.5 * K, length=101)
h_grid <- x_grid
```

With parameters `r p`.

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

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

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

Which estimates a Ricker model with $r =$ `r alt$p[1]`, $K =$ `r alt$p[2]`, and the Allen allee model with $r =$ `r est$p[1]`, $K =$ `r est$p[2]` and $C =$ `r est$p[3]`.

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

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

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

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


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

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

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

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

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

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

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

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

## ----sdp-pars, dependson="stateeq"---------------------------------------
sigma_g <- 0.01
sigma_m <- 0.0
z_g <- function() rlnorm(1, 0, sigma_g)
z_m <- function() 1
x_grid <- seq(0, 1.5 * K, length=50)
h_grid <- x_grid
profit <- function(x,h) pmin(x, h)
delta <- 0.01
OptTime <- 50  # stationarity with unstable models is tricky thing
reward <- 0
xT <- 0
Xo <-  allee+.5# observations start from
x0 <- 10 # simulation under policy starts from
Tobs <- 40
MaxT = 1000 # timeout for value iteration convergence

## ----obs, dependson="sdp-pars", fig.keep='none'--------------------------
  set.seed(1234)
  #harvest <- sort(rep(seq(0, .5, length=7), 5))
  x <- numeric(Tobs)
  x[1] <- Xo
  nz <- 1
  for(t in 1:(Tobs-1))
    x[t+1] = z_g() * f(x[t], h=0, p=p)
  obs <- data.frame(x = c(rep(0,nz),
                          pmax(rep(0,Tobs-1), x[1:(Tobs-1)])),
                    y = c(rep(0,nz),
                          x[2:Tobs]))
raw_plot <- ggplot(data.frame(time = 1:Tobs, x=x), aes(time,x)) + geom_line()

## ----mle, dependson="obs"------------------------------------------------
set.seed(12345)
estf <- function(p){
    mu <- f(obs$x,0,p)
    -sum(dlnorm(obs$y, log(mu), p[4]), log=TRUE)
}
par <- c(p[1]*rlnorm(1,0,.1),
         p[2]*rlnorm(1,0,.1),
         p[3]*rlnorm(1,0, .1),
         sigma_g * rlnorm(1,0,.1))
o <- optim(par, estf, method="L", lower=c(1e-5,1e-5,1e-5,1e-5))
f_alt <- f
p_alt <- c(as.numeric(o$par[1]), as.numeric(o$par[2]), as.numeric(o$par[3]))
sigma_g_alt <- as.numeric(o$par[4])

est <- list(f = f_alt, p = p_alt, sigma_g = sigma_g_alt, mloglik=o$value)

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

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

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

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

## ----gp-output, dependson="gp", results=toggle---------------------------
tgp_dat <-
    data.frame(  x = x_grid,
                 y = gp_dat$E_Ef,
                 ymin = gp_dat$E_Ef - 2 * sqrt(gp_dat$E_Vf),
                 ymax = gp_dat$E_Ef + 2 * sqrt(gp_dat$E_Vf) )

## ----jags-setup----------------------------------------------------------
y <- x
N <- length(x);
jags.data <- list("N","y")
n.chains <- 6
n.iter <- 1e6
n.burnin <- floor(10000)
n.thin <- max(1, floor(n.chains * (n.iter - n.burnin)/1000))
n.update <- 10

## ----common-priors-------------------------------------------------------
stdQ_prior_p <- c(1e-6, 100)
stdR_prior_p <- c(1e-6, .1)
stdQ_prior  <- function(x) dunif(x, stdQ_prior_p[1], stdQ_prior_p[2])
stdR_prior  <- function(x) dunif(x, stdR_prior_p[1], stdR_prior_p[2])

## ----allen-model---------------------------------------------------------
K_prior_p <- c(0.01, 20.0)
r0_prior_p <- c(0.01, 6.0)
theta_prior_p <- c(0.01, 20.0)

bugs.model <-
paste(sprintf(
"model{
  K     ~ dunif(%s, %s)
  r0    ~ dunif(%s, %s)
  theta ~ dunif(%s, %s)
  stdQ ~ dunif(%s, %s)",
  K_prior_p[1], K_prior_p[2],
  r0_prior_p[1], r0_prior_p[2],
  theta_prior_p[1], theta_prior_p[2],
  stdQ_prior_p[1], stdQ_prior_p[2]),

  "
  iQ <- 1 / (stdQ * stdQ);
  y[1] ~ dunif(0, 10)
  for(t in 1:(N-1)){
    mu[t] <- log(y[t]) + r0 * (1 - y[t]/K)* (y[t] - theta) / K
    y[t+1] ~ dlnorm(mu[t], iQ)
  }
}")
writeLines(bugs.model, "allen_process.bugs")

## ----allen-priors, dependson="common-priors"-----------------------------
K_prior     <- function(x) dunif(x, K_prior_p[1], K_prior_p[2])
r0_prior <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2])
theta_prior <- function(x) dunif(x, theta_prior_p[1], theta_prior_p[2])
par_priors  <- list(K = K_prior, deviance = function(x) 0 * x,
                    r0 = r0_prior, theta = theta_prior,
                    stdQ = stdQ_prior)

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

## ----allen-mcmc, dependson=c("allen-model", "allen-pars", "jags-setup"), results=toggle----
jags.params=c("K","r0","theta","stdQ") # be sensible about the order here
jags.inits <- function(){
  list("K"= 10 * rlnorm(1,0, 0.1),
       "r0"= 1 * rlnorm(1,0, 0.1) ,
       "theta"=   5 * rlnorm(1,0, 0.1) ,
       "stdQ"= abs( 0.1 * rlnorm(1,0, 0.1)),
       .RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}

set.seed(1234)
# parallel refuses to take variables as arguments (e.g. n.iter = 1e5 works, but n.iter = n doesn't)
allen_jags <- do.call(jags, list(data=jags.data, inits=jags.inits,
                                      jags.params, n.chains=n.chains,
                                      n.iter=n.iter, n.thin=n.thin,
                                      n.burnin=n.burnin,
                                      model.file="allen_process.bugs"))

# Run again iteratively if we haven't met the Gelman-Rubin convergence criterion
recompile(allen_jags) # required for parallel
allen_jags <- do.call(autojags,
											list(object=allen_jags, n.update=n.update,
                           n.iter=n.iter, n.thin = n.thin))

## ----allen-traces, dependson="allen-mcmc"--------------------------------
tmp <- lapply(as.mcmc(allen_jags), as.matrix) # strip classes the hard way...
allen_posteriors <- melt(tmp, id = colnames(tmp[[1]]))
names(allen_posteriors) = c("index", "variable", "value", "chain")
plot_allen_traces <- ggplot(allen_posteriors) + geom_line(aes(index, value)) +
  facet_wrap(~ variable, scale="free", ncol=1)

## ----allen-posteriors, dependson=c("allen-traces", "allen-priors")-------
allen_priors <- ddply(allen_posteriors, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})
plot_allen_posteriors <- ggplot(allen_posteriors, aes(value)) +
  stat_density(geom="path", position="identity", alpha=0.7) +
#  geom_line(data=allen_priors, aes(x=value, y=density), col="red") +
  facet_wrap(~ variable, scale="free", ncol=3)

## ----allen-output, dependson=c("posterior-mode", "allen-traces"), results=toggle----
A <- allen_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
pardist <- acast(A, index ~ variable)
bayes_coef <- apply(pardist,2, posterior.mode)
bayes_pars <- unname(c(bayes_coef["r0"], bayes_coef["K"], bayes_coef["theta"])) # parameters formatted for f
allen_f <- function(x,h,p) unname(RickerAllee(x,h, unname(p[c("r0", "K", "theta")])))
allen_means <- sapply(x_grid, f, 0, bayes_pars)
bayes_pars
head(pardist)

## ----ricker-model--------------------------------------------------------
K_prior_p <- c(0.01, 40.0)
r0_prior_p <- c(0.01, 20.0)
bugs.model <-
paste(sprintf(
"model{
  K    ~ dunif(%s, %s)
  r0    ~ dunif(%s, %s)
  stdQ ~ dunif(%s, %s)",
  K_prior_p[1], K_prior_p[2],
  r0_prior_p[1], r0_prior_p[2],
  stdQ_prior_p[1], stdQ_prior_p[2]),
  "
  iQ <- 1 / (stdQ * stdQ);
  y[1] ~ dunif(0, 10)
  for(t in 1:(N-1)){
    mu[t] <- log(y[t]) + r0 * (1 - y[t]/K)
    y[t+1] ~ dlnorm(mu[t], iQ)
  }
}")
writeLines(bugs.model, "ricker_process.bugs")

## ----ricker-priors, dependson="common-priors"----------------------------
K_prior     <- function(x) dunif(x, K_prior_p[1], K_prior_p[2])
r0_prior <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2])
par_priors <- list(K = K_prior, deviance = function(x) 0 * x,
                   r0 = r0_prior, stdQ = stdQ_prior)

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

## ----ricker-mcmc, dependson="ricker-model", results=toggle---------------
jags.params=c("K","r0", "stdQ")
jags.inits <- function(){
  list("K"= 10 * rlnorm(1,0,.5),
       "r0"= rlnorm(1,0,.5),
       "stdQ"=sqrt(0.05) * rlnorm(1,0,.5),
       .RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
set.seed(12345)
ricker_jags <- do.call(jags,
                       list(data=jags.data, inits=jags.inits,
                            jags.params, n.chains=n.chains,
                            n.iter=n.iter, n.thin=n.thin, n.burnin=n.burnin,
                            model.file="ricker_process.bugs"))
recompile(ricker_jags)
ricker_jags <- do.call(autojags,
                       list(object=ricker_jags, n.update=n.update,
														n.iter=n.iter, n.thin = n.thin,
														progress.bar="none"))

## ----ricker-traces, dependson="ricker-mcmc"------------------------------
tmp <- lapply(as.mcmc(ricker_jags), as.matrix) # strip classes the hard way...
ricker_posteriors <- melt(tmp, id = colnames(tmp[[1]]))
names(ricker_posteriors) = c("index", "variable", "value", "chain")
plot_ricker_traces <- ggplot(ricker_posteriors) + geom_line(aes(index, value)) +
  facet_wrap(~ variable, scale="free", ncol=1)

## ----ricker-posteriors, dependson=c("ricker-traces", "ricker-priors")----
ricker_priors <- ddply(ricker_posteriors, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})
# plot posterior distributions
plot_ricker_posteriors <- ggplot(ricker_posteriors, aes(value)) +
  stat_density(geom="path", position="identity", alpha=0.7) +
#  geom_line(data=ricker_priors, aes(x=value, y=density), col="red") +  # don't plot priors
  facet_wrap(~ variable, scale="free", ncol=2)

## ----ricker-output, dependson=c("posterior-mode", "ricker-traces"), results=toggle----
A <- ricker_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
ricker_pardist <- acast(A, index ~ variable)
bayes_coef <- apply(ricker_pardist,2, posterior.mode)
ricker_bayes_pars <- unname(c(bayes_coef["r0"], bayes_coef["K"]))
ricker_f <- function(x,h,p){
  sapply(x, function(x){
    x <- pmax(0, x-h)
    pmax(0, x * exp(p["r0"] * (1 - x / p["K"] )) )
  })
}
ricker_means <- sapply(x_grid, Ricker, 0, ricker_bayes_pars[c(1,2)])
head(ricker_pardist)
ricker_bayes_pars

## ----myers-model---------------------------------------------------------
r0_prior_p <- c(.0001, 10.0)
theta_prior_p <- c(.0001, 10.0)
K_prior_p <- c(.0001, 40.0)
bugs.model <-
paste(sprintf(
"model{
  r0    ~ dunif(%s, %s)
  theta    ~ dunif(%s, %s)
  K    ~ dunif(%s, %s)
  stdQ ~ dunif(%s, %s)",
  r0_prior_p[1], r0_prior_p[2],
  theta_prior_p[1], theta_prior_p[2],
  K_prior_p[1], K_prior_p[2],
  stdQ_prior_p[1], stdQ_prior_p[2]),

  "
  iQ <- 1 / (stdQ * stdQ);

  y[1] ~ dunif(0, 10)
  for(t in 1:(N-1)){
    mu[t] <- log(r0)  + theta * log(y[t]) - log(1 + pow(abs(y[t]), theta) / K)
    y[t+1] ~ dlnorm(mu[t], iQ)
  }
}")
writeLines(bugs.model, "myers_process.bugs")

## ----myers-priors--------------------------------------------------------
K_prior     <- function(x) dunif(x, K_prior_p[1], K_prior_p[2])
r_prior     <- function(x) dunif(x, r0_prior_p[1], r0_prior_p[2])
theta_prior <- function(x) dunif(x, theta_prior_p[1], theta_prior_p[2])
par_priors <- list( deviance = function(x) 0 * x, K = K_prior,
                    r0 = r_prior, theta = theta_prior,
                    stdQ = stdQ_prior)

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

## ----myers-mcmc, dependson="myers-model", results=toggle-----------------
jags.params=c("r0", "theta", "K", "stdQ")
jags.inits <- function(){
  list("r0"= 1 * rlnorm(1,0,.1),
       "K"=    10 * rlnorm(1,0,.1),
       "theta" = 1 * rlnorm(1,0,.1),
       "stdQ"= sqrt(0.2) * rlnorm(1,0,.1),
       .RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
set.seed(12345)
myers_jags <- do.call(jags,
                      list(data=jags.data, inits=jags.inits,
													 jags.params, n.chains=n.chains,
													 n.iter=n.iter, n.thin=n.thin,
                           n.burnin=n.burnin,
                           model.file="myers_process.bugs"))
recompile(myers_jags)
myers_jags <- do.call(autojags,
                      list(myers_jags, n.update=n.update,
                           n.iter=n.iter, n.thin = n.thin,
                           progress.bar="none"))

## ----myers-traces, dependson="myers-mcmc"--------------------------------
tmp <- lapply(as.mcmc(myers_jags), as.matrix) # strip classes
myers_posteriors <- melt(tmp, id = colnames(tmp[[1]]))
names(myers_posteriors) = c("index", "variable", "value", "chain")
plot_myers_traces <- ggplot(myers_posteriors) +
  geom_line(aes(index, value)) + # priors, need to fix order though
  facet_wrap(~ variable, scale="free", ncol=1)

## ----myers-posteriors, dependson="myers-traces"--------------------------
par_prior_curves <- ddply(myers_posteriors, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})
plot_myers_posteriors <- ggplot(myers_posteriors, aes(value)) +
  stat_density(geom="path", position="identity", alpha=0.7) +
#  geom_line(data=par_prior_curves, aes(x=value, y=density), col="red") +  # Whoops, these are misaligned. see table instead
  facet_wrap(~ variable, scale="free", ncol=3)

## ----myers-output, dependson=c("posterior-mode", "myers-traces"), results=toggle----
A <- myers_posteriors
A$index <- A$index + A$chain * max(A$index) # Combine samples across chains by renumbering index
myers_pardist <- acast(A, index ~ variable)
bayes_coef <- apply(myers_pardist,2, posterior.mode) # much better estimates
myers_bayes_pars <- unname(c(bayes_coef["r0"], bayes_coef["theta"], bayes_coef["K"]))
myers_means <- sapply(x_grid, Myer_harvest, 0, myers_bayes_pars)
myers_f <- function(x,h,p) Myer_harvest(x, h, p[c("r0", "theta", "K")])
head(myers_pardist)
myers_bayes_pars

## ----assemble-models, dependson=c("myers-output", "ricker-output", "allen-output", "gp-output", "mle-output")----
models <- data.frame(x=x_grid,
										 GP=tgp_dat$y,
										 True=true_means,
                     MLE=est_means,
										 Ricker=ricker_means,
                     Allen = allen_means,
                     Myers = myers_means)
models <- melt(models, id="x")
# some labels
names(models) <- c("x", "method", "value")
# labels for the colorkey too
model_names = c("GP", "True", "MLE", "Ricker", "Allen", "Myers")
colorkey=cbPalette
names(colorkey) = model_names
step_ahead <- function(x, f, p){
  h = 0
  x_predict <- sapply(x, f, h, p)
  n <- length(x_predict) - 1
  y <- c(x[1], x_predict[1:n])
  y
}

## ----gp-opt, dependson="gp-output"---------------------------------------
matrices_gp <- gp_transition_matrix(gp_dat$Ef_posterior, gp_dat$Vf_posterior, x_grid, h_grid)
opt_gp <- value_iteration(matrices_gp, x_grid, h_grid, MaxT, xT, profit, delta, reward)

## ----mle-opt, dependson="mle-output"-------------------------------------
matrices_true <- f_transition_matrix(f, p, x_grid, h_grid, sigma_g)
opt_true <- value_iteration(matrices_true, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
matrices_estimated <- f_transition_matrix(est$f, est$p, x_grid, h_grid, est$sigma_g)
opt_estimated <- value_iteration(matrices_estimated, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)

## ----allen-opt, dependson="allen-output"---------------------------------
matrices_allen <- parameter_uncertainty_SDP(allen_f, x_grid, h_grid, pardist, 4)
opt_allen <- value_iteration(matrices_allen, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)

## ----ricker-opt, dependson="ricker-output"-------------------------------
matrices_ricker <- parameter_uncertainty_SDP(ricker_f, x_grid, h_grid, as.matrix(ricker_pardist), 3)
opt_ricker <- value_iteration(matrices_ricker, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)

## ----myers-opt, dependson="myers-output"---------------------------------
matrices_myers <- parameter_uncertainty_SDP(myers_f, x_grid, h_grid, as.matrix(myers_pardist), 4)
myers_alt <- value_iteration(matrices_myers, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)

## ----assemble-opt, dependson=c("gp-opt", "mle-opt", "allen-opt", "ricker-opt", "myers-opt")----
OPT = data.frame(GP = opt_gp$D, True = opt_true$D, MLE = opt_estimated$D, Ricker = opt_ricker$D, Allen = opt_allen$D, Myers = myers_alt$D)
colorkey=cbPalette
names(colorkey) = names(OPT)

## ----sims, dependson="assemble-opt"--------------------------------------
sims <- lapply(OPT, function(D){
  set.seed(1)
  lapply(1:100, function(i)
    ForwardSimulate(f, p, x_grid, h_grid, x0, D, z_g, profit=profit, OptTime=OptTime)
  )
})
# turn the list into a data.frame
keep.as.columns <- names(sims[[1]][[1]])
dat <- melt(sims, id=keep.as.columns)
sims_data <- data.table(dat)
setnames(sims_data, c("L1", "L2"), c("method", "reps"))
# Legend in original ordering please, not alphabetical:
sims_data$method = factor(sims_data$method, ordered=TRUE, levels=names(OPT))

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

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

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

## ----deviances, dependson=c("dic_calc"), include=FALSE, echo=FALSE-------
allen_deviance  <- - posterior.mode(pardist[,'deviance'])
ricker_deviance <- - posterior.mode(ricker_pardist[,'deviance'])
myers_deviance  <- - posterior.mode(myers_pardist[,'deviance'])
true_deviance   <- 2*estf(c(p, sigma_g))
mle_deviance    <- 2*estf(c(est$p, est$sigma_g))
aictable <- data.frame(Allen = allen_deviance + 2*(1+length(bayes_pars)),  # +1 for noise parameter
                       Ricker = ricker_deviance + 2*(1+length(ricker_bayes_pars)),
                       Myers = myers_deviance + 2*(1+length(myers_bayes_pars)),
                       row.names = c("AIC"))
bictable <- data.frame(Allen = allen_deviance + log(length(x))*(1+length(bayes_pars)),
                       Ricker = ricker_deviance + log(length(x))*(1+length(ricker_bayes_pars)),
                       Myers = myers_deviance + log(length(x))*(1+length(myers_bayes_pars)),
                       row.names = c("BIC"))
xtable::xtable(rbind(dictable, aictable, bictable))

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

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

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

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

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

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

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

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