text:
\documentclass[author-year, 12pt,review]{elsarticle} %review=doublespace preprint=single 5p=2 column
%%% Begin My package additions %%%%%%%%%%%%%%%%%%%
\usepackage[hyphens]{url}
\usepackage{lineno} % add
  \linenumbers % turns line numbering on
\bibliographystyle{elsarticle-harv}
\biboptions{sort&compress} % For natbib
\usepackage{graphicx}
\usepackage{booktabs} % book-quality tables
%% Redefines the elsarticle footer
\makeatletter
\def\ps@pprintTitle{%
 \let\@oddhead\@empty
 \let\@evenhead\@empty
 \def\@oddfoot{\it \hfill\today}%
 \let\@evenfoot\@oddfoot}
\makeatother

% A modified page layout
\textwidth 6.75in
\oddsidemargin -0.15in
\evensidemargin -0.15in
\textheight 9in
\topmargin -0.5in
%%%%%%%%%%%%%%%% end my additions to header

\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
  \usepackage[utf8]{inputenc}
\else % if luatex or xelatex
  \usepackage{fontspec}
  \ifxetex
    \usepackage{xltxtra,xunicode}
  \fi
  \defaultfontfeatures{Mapping=tex-text,Scale=MatchLowercase}
  \newcommand{\euro}{€}
\fi
% use microtype if available
\IfFileExists{microtype.sty}{\usepackage{microtype}}{}
\usepackage{longtable}
\usepackage{graphicx}
% We will generate all images so they have a width \maxwidth. This means
% that they will get their normal width if they fit onto the page, but
% are scaled down if they would overflow the margins.
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth
\else\Gin@nat@width\fi}
\makeatother
\let\Oldincludegraphics\includegraphics
\renewcommand{\includegraphics}[1]{\Oldincludegraphics[width=\maxwidth]{#1}}
\ifxetex
  \usepackage[setpagesize=false, % page size defined by xetex
              unicode=false, % unicode breaks when used with xetex
              xetex]{hyperref}
\else
  \usepackage[unicode=true]{hyperref}
\fi
\hypersetup{breaklinks=true,
            bookmarks=true,
            pdfauthor={},
            pdftitle={Avoiding tipping points in fisheries management through Gaussian Process Dynamic Programming},
            colorlinks=true,
            urlcolor=blue,
            linkcolor=magenta,
            pdfborder={0 0 0}}
\urlstyle{same}  % don't use monospace font for urls
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
\setlength{\emergencystretch}{3em}  % prevent overfull lines
\setcounter{secnumdepth}{0}
% Pandoc toggle for numbering sections (defaults to be off)
\setcounter{secnumdepth}{0}
% Pandoc header

\begin{document}
\begin{frontmatter}

  \title{Avoiding tipping points in fisheries management through Gaussian Process
Dynamic Programming}
    \author[cstar]{Carl Boettiger\corref{c1}}
   \ead{cboettig(at)gmail.com}
   \cortext[c1]{Corresponding author}
    \author[cstar]{Marc Mangel}

    \author[noaa]{Stephan Munch}

      \address[cstar]{Center for Stock Assessment Research, Department of Applied Math and
Statistics, University of California, Mail Stop SOE-2, Santa Cruz, CA
95064, USA}
    \address[marc]{Center for Stock Assessment Research, Department of Applied Math and
Statistics, University of California, Mail Stop SOE-2, Santa Cruz, CA
95064, USA and Department of Biology, University of Bergen, Bergen,
Norway 9020}
    \address[noaa]{Southwest Fisheries Science Center, National Oceanic and Atmospheric
Administration, 110 Shaffer Road, Santa Cruz, CA 95060, USA}

  \begin{abstract}
  Model uncertainty and limited data are fundamental challenges to robust
  management of human intervention in a natural system. These challenges
  are acutely highlighted by concerns that many ecological systems may
  contain tipping points, such as Allee population sizes. Before a
  collapse, we do not know where the tipping points lie, if they exist at
  all. Hence, we know neither a complete model of the system dynamics nor
  do we have access to data in some large region of state-space where such
  a tipping point might exist. We illustrate how a Bayesian Non-Parametric
  (BNP) approach using a Gaussian Process (GP) prior provides a flexible
  representation of this inherent uncertainty. We embed GPs in a
  Stochastic Dynamic Programming (SDP) framework in order to make robust
  management predictions with both model uncertainty and limited data. We
  use simulations to evaluate this approach as compared with the standard
  approach of using model selection to choose from a set of candidate
  models. We find that model selection erroneously favors models without
  tipping points -- leading to harvest policies that guarantee extinction.
  The GPDP performs nearly as well as the true model and significantly
  outperforms standard approaches. We illustrate this using examples of
  simulated single-species dynamics, where the standard model selection
  approach should be most effective, and find that it still fails to
  account for uncertainty appropriately and leads to population crashes,
  while management based on the GPDP does not, since it does not
  underestimate the uncertainty outside of the observed data.
  \end{abstract}
   \begin{keyword} Bayesian \sep Structural Uncertainty \sep Nonparametric \sep Optimal Control \sep Decision Theory \sep Gaussian Processes \sep Fisheries Management \sep \end{keyword}
 \end{frontmatter}

\section{Introduction}\label{introduction}

Decision making under uncertainty is a ubiquitous challenge in the
management of human intervention in natural resources and conservation.
Decision-theoretic approaches provide a framework to determine the best
sequence of actions in face of uncertainty, but only when that
uncertainty can be meaningfully quantified (Fischer et al. 2009). Over
the last four decades (beginning with Clark (1976), Clark (2009) and
Walters and Hilborn (1978)) dynamic optimization methods, particularly
Stochastic Dynamic Programming (SDP), have become increasingly important
as a means of understanding how to manage human intervention into
natural systems. Simultaneously, there has been increasing recognition
of the importance of multiple steady states or `tipping points'
(Scheffer et al. 2001, 2009, Polasky et al. 2011) in ecological systems.

We develop a novel approach to address these concerns in the context of
fisheries; although the challenges and methods are germane to other
problems of conservation or natural resource exploitation. Economic
value and ecological concern have made marine fisheries the crucible for
much of the founding work on management under uncertainty (Gordon 1954,
Clark 1976, 2009, May et al. 1979, Reed 1979, Ludwig and Walters 1982).

Even if we know the proper deterministic description of the biological
system, there is intrinsic stochasticity in biological dynamics,
measurements, and implementation of policy (\emph{e.g.} Reed 1979, Clark
and Kirkwood 1986, Roughgarden and Smith 1996, Sethi et al. 2005). We
may also lack knowledge about the parameters of the biological dynamics
(parametric uncertainty, \emph{e.g.} Ludwig and Walters 1982, Hilborn
and Mangel 1997, McAllister 1998, Schapaugh and Tyre 2013), or even not
know which model is proper description of the system (structural
uncertainty, \emph{e.g.} Williams 2001, Cressie et al. 2009,
Athanassoglou and Xepapadeas 2012). Of these, the latter is generally
the hardest to quantify. Typical approaches confront the data with a
collection of models, assuming that the true dynamics (or reasonable
approximation) is among the collection and then use model choice or
model averaging to arrive at a conclusion (Williams 2001, Cressie et al.
2009, Athanassoglou and Xepapadeas 2012). Even setting aside other
concerns (see Cressie et al. (2009)), these approaches are unable to
describe uncertainty outside the observed data range.

Structural uncertainty is particularly insidious when we try to predict
outside of the range of observed data (Mangel et al. 2001) because we
are extrapolating into unknown regions. In management applications, this
extrapolation uncertainty is particularly important since (a) management
involves considering actions that may move the system outside the range
of observed behavior, and (b) the decision tools (optimal control
theory, SDP) rely on both reasonable estimates of the expected outcomes
and on the weights given to those outcomes (\emph{e.g.} Weitzman 2013).
Thus characterizing uncertainty is as important as characterizing the
expected outcome.

Tipping points in ecological dynamics (Scheffer et al. 2001, Polasky et
al. 2011) highlight this problem because precise models are not
available and data are limited such as around high stock levels or an
otherwise desirable state. With perfect information, one would know just
how far a system could be pushed before crossing the tipping point, and
management would be simple. But we face imperfect models and limited
data and, with tipping points,\\even small errors can have very large
consequences, as we shall illustrate later. Because intervention may be
too late once a tipping point has been crossed (but see Hughes et al.
(2013)), management is often concerned with avoiding tipping points
before any data about them are available.

The dual concerns of model uncertainty and incomplete data create a
substantial challenge to existing decision-theoretic approaches
(Brozović and Schlenker 2011). We illustrate how Stochastic Dynamic
Programming (SDP) (Mangel and Clark 1988, Marescot et al. 2013) can be
implemented using a Bayesian Non-Parametric (BNP) model of population
dynamics (Munch et al. 2005a). The BNP method has two distinct
advantages. First, using a BNP model sidesteps the need for an accurate
model-based description of the system dynamics. Second, a BNP model
better reflects uncertainty when extrapolating beyond the observed data.
This is crucial to providing robust decision-making when the correct
model is not known (as is almost always true). {[}We use \emph{robust}
to characterize approaches that provide nearly optimal solutions without
being sensitive to the choice of the (unknown) underlying model.{]}

This paper is the first ecological application of the SDP without an
\emph{a priori} model of the underlying dynamics. Unlike parametric
approaches that can only reflect uncertainty in parameter estimates, the
BNP method provides a broader representation of uncertainty, including
uncertainty beyond the observed data. We will show that Gaussian Process
Dynamic Programming (GPDP) allows us to find robust management solutions
in face of limited data without knowing the correct model structure.

For comparisons, we consider the performance of management based on GPDP
against management policies derived under several alternative parametric
models (Reed 1979, Ludwig and Walters 1982, Mangel and Clark 1988).
Rather than compare models in terms of best fit to data, we compare
model performance in the concrete terms of the decision-maker's
objectives.

\section{Approach and Methods}\label{approach-and-methods}

We first describe the requirements of dynamic optimization for the
management of human intervention in natural resource systems. After that
we describe three parametric models for population dynamics and the
Gaussian Process (GP)\footnote{We abbreviate Gaussian Process as GP,
  which refers to the statistical model we use to approximate the
  population dynamics, and we use the term Gaussian Process Dynamic
  Programming (GPDP), to refer to the \emph{use} of a GP as the
  underlying process model when solving a Dynamic Programming equation.
  Hence we will refer to the models as: GP, Ricker, Allen, etc, and the
  novel method we put forward here as GPDP.} description of population
dynamics. All computer code used here has been embedded in the
manuscript sources (see Xie (2013)), and an implementation of the GPDP
approach is provided as an accompanying R package. Source code, R
package and the CSV data files corresponding each figure are archived in
the supplement (Boettiger et al. 2014).

\subsubsection{Requirements of dynamic
optimization}\label{requirements-of-dynamic-optimization}

Dynamic optimization equires characterizing the dynamics of a state
variable (or variables), a control action, and a value function. For
simplicity, we consider only a single state variable. This is a
best-case scenario for the parametric models because we simulate
underlying dynamics from one of the three parametric models , whereas in
the natural world we never know the ``true'' model. In addition, by
choosing one-dimensional models with just a few parameters, we limit the
chance that poor performance will be due to inability to estimate
parameters accurately, something that becomes a more severe problem for
higher-dimensional parametric models. Finally, the parametric models we
consider are commonly used in modeling stock-recruitment dynamics or to
model sudden transitions between alternative stable states.

We let \(X(t)\) denote the size (numbers or biomass) of the focal
population at time \(t\) and assume that in the absence of take its
dynamics are:

\begin{equation}
X(t+1) = Z(t) f(X(t), \mathbf{p}) \label{eq1}
\end{equation}

Where Z(t) is log-normally distributed process stochasticity (Reed 1979)
and \(\mathbf{p}\) is a vector of parameters to be estimated from the
data. We describe the three choices for \(f(X(t),\mathbf{p})\) in the
next section.\\The control action is a harvest or take, \(h(t)\),
measured in the same units as \(X\), at time \(t\). Thus, in the
presence of take, the population size on the right hand side of Eqn 1 is
replaced by \(S(t) = X(t) - h(t)\).

To construct the value function, we consider a return when
\(X(t) = x(t)\) and harvest \(h(t) = h\) denoted as the reward,
\(R(x(t), h)\). For example, if the return is the harvest at time \(t\),
then \(R(x(t), h(t)) = \min(x(t), h(t))\). We assume that future
harvests are discounted relative to current ones at a constant rate of
discount \(\delta\) and ask for the harvest policy that maximizes total
discounted harvest between the current time \(t\) and a final time
\(T\). That is, we seek to maximize over choices of harvest
\(E [ \sum_{t = 0}^{T} R(X(t), h(t), t) \delta^t]\), where the state
dynamics are given by Eqn 1 and \(E\) denotes the expectation over
future population states.

In order to find that policy, we introduce the value function
\(V(x(t), t)\) representing the total discounted catch from time \(t\)
onwards given that \(X(t) = x(t)\). This value function satisfies an
equation of SDP (Mangel and Clark 1988, Clark and Mangel 2000, Clark
2009 , Mangel 2014),

\begin{equation}
V(x(t), t) = \max_{h(t)} \lbrace R(h(t), x(t)) + \delta \cdot \mathbf{\mathrm{E}}_{X(t+1)} \left[ V(X(t+1), t+1) | x(t), h(t) \right] \rbrace
\end{equation}

where expectation is taken over all possible values of the next state,
\(X(t+1)\), and maximized over all possible choices of harvest,
\(h(t)\). That is, at time \(t\), when population size is \(x(t)\) and
harvest \(h(t)\) is applied, the immediate return is \(R(h(t), x(t))\).
When the sole source of uncertainty is the process stochasticity term,
\(Z\), and thus the expectation in Eqn 2 is equivalent to taking
expectations over \(Z(t)\). That is

\begin{equation}
\mathbf{\mathrm{E}}_{X(t+1)} \left[ V( X(t+1),t+1) | x(t), h(t) \right] = \mathbf{\mathrm{E}}_{Z(t)} \left[V( Z(t) f(x(t) - h(t))|\mathbf{p}), t+1 | x(t), h(t) \right]
\end{equation}

where the population size after the take is \(x(t) - h(t)\), which is
then translated into \(X(t+1)\) by Eqn 1 (that is, we replace \(X(t+1)\)
by \(Z(t) f(x(t)-h(t)|\mathbf{p})\)).

When the parameters governing the dynamics are also uncertain, we take
the expectation over the posterior distribution for the parameters:

\begin{equation}
\mathbf{\mathrm{E}}_{X(t+1)} \left[ V(X(t+1), t+1) | x(t), h(t) \right] = \mathbf{\mathrm{E}}_{\mathbf{p}|\mathrm{data}} \{ \mathbf{\mathrm{E}}_{Z(t) | \mathbf{p}, \mathrm{data}} \left[ V(Z(t) f(x(t) - h(t)|\mathbf{p}), t)  \right] \}
\end{equation}

When the underlying population dynamics are unknown (the case of
structural uncertainty), the function \(f\) itself is uncertain and the
expectation for the next state includes uncertainty in \(f\) as well.
That is

\begin{equation}
\mathbf{\mathrm{E}}_{X(t+1)} \left[ V(X(t+1), t+1) | x(t), h(t) \right] = \mathbf{\mathrm{E}}_{\mathbf{p}|\mathrm{data}} \{ \mathbf{\mathrm{E}}_{f, Z(t) | \mathbf{p}, \mathrm{data}} \left[ V( Z(t) f(x(t) - h(t)| \mathbf{p}), t) \right] \}
\end{equation}

We consider the finite time problem with \(T=\) 1000, which we solve
using the standard value iteration algorithm (see Mangel and Clark 1988,
Clark and Mangel 2000).

\subsubsection{Parametric Models}\label{parametric-models}

We consider three candidate parametric models for the population
dynamics: The Ricker model, the Allen model (Allen and Tanner 2005), and
the Myers model (Myers et al. 1995), Eqns \eqref{ricker}-\eqref{myers}.
In all three, we let \(K\) denote the carrying capacity and \(r\) the
maximum per capita growth rate. The Ricker model has two parameters and
the right hand side of Eqn 1 is

\begin{equation}
f(S(t)|r,K) = S(t) e^{r \left(1 - \frac{S(t)}{K} \right) } \label{ricker}
\end{equation}

The Allen model has three parameters

\begin{equation}
f(S(t)|r, K, X_c) = S(t) e^{r \left(1 - \frac{S(t)}{K}\right)\left(S(t) - X_c\right)} \label{allen}
\end{equation}

where \(X_c\) denotes the location of the unstable steady state (i.e.,
the tipping point).

The Myers model also has three parameters

\begin{equation}
f(S(t) | r, K, \theta)  = \frac{r S(t)^{\theta}}{1 + \frac{S(t)^\theta}{K}} \label{myers}
\end{equation}

where \(\theta = 1\) corresponds to Beverton-Holt dynamics and
\(\theta\) \textgreater{} 2 leads to Allee effects and multiple stable
states.

The Ricker model does not lead to multiple steady states. The Allen
model resembles the Ricker dynamics with an added Allee effect parameter
(Courchamp et al. 2008), below which the population cannot persist. The
Myers model also has three parameters and contains an Allee threshold,
but unlike the Ricker model saturates at high population size. The
multiplicative log-normal stochasticity in Eqn 1 introduces one
additional parameter \(\sigma\) that must be estimated.

Because of our interest in management performance in the presence of
tipping points, all of our simulations are based on the Allen model. The
Allen model is thus the state of nature and is expected to provide the
best-case scenario. The Ricker model is a reasonable approximation of
these dynamics far from the Allee threshold (but lacks threshold
dynamics), while the Myers model shares the essential feature of a
threshold but differs in structure from the Allen model. Throughout, we
refer to the ``True'' model when the underlying parameters \emph{are
known without error}, and refer to the ``Allen'' model when these
parameters have been estimated from the sample data.

We consider a period of 40 in which data are obtained to estimate the
parameters or the GP. This is long enough that the estimates do not
depend on the particular realization, and longer times are not likely to
provide substantial improvement. Each of the models is fit to the same
data (Figure 1).

We inferred posterior distributions for the parameters of each model by
Gibbs sampling (Gelman et al. (2003) implemented in R (R Core Team 2013)
using \texttt{jags}, (Su and Masanao Yajima 2013)). We choose uniform
priors for all parameters of the parametric models (See appendix Tables
S1-S3; R code provided). We show one-step-ahead predictions of these
model fits in Figure 1. We tested each chain for Gelman-Rubin
convergence and results were robust to longer runs. For each simulation
we also applied several commonly used model selection criteria (AIC,
BIC, DIC, see Burnham and Anderson (2002)) to identify the best fitting
model.

Additionally, we compute the maximum likelihood estimate (MLE, as we
will refer to this model in the figures) of the parameters for the
(structurally correct) Allen model. Comparing this to using the
posterior distribution of parameters inferred from MCMC for the same
model gives some indication of the importance of this uncertainty in the
dynamic programming.

\subsubsection{The Gaussian Process
model}\label{the-gaussian-process-model}

The core difference for our purpose between the estimated GP and the
estimated parametric models is that the estimated GP model is defined
explicitly in reference to the observed data. As a result, uncertainty
arises in the GP model not only from uncertainty in the parameters, but
is also increases in regions farther from the observed states, such as
low population sizes in the example illustrated here. The estimated
parametric models, by contrast, are completely specified by the
parameters.

The use of GPs to characterize dynamical systems is relatively new
(Kocijan et al. 2005), and was first introduced in the context
ecological modeling and fisheries management in Munch et al. (2005b). GP
models have subsequently been used to test for the presence of Allee
effects (Sugeno and Munch 2013a), estimate the maximum reproductive rate
(Sugeno and Munch 2013b), determine temporal variation in food
availability (Sigourney et al. 2012), and provide a basis for
identifying model-misspecification (Thorson et al. 2014). An accessible
and thorough introduction to the formulation and use of GPs can be found
in Rasmussen and Williams (2006).

A GP is a stochastic process for which any realization consisting of n
points follows a multivariate normal distribution of dimension \(n\). To
characterize the GP we need a mean function and a covariance function.
We proceed as follows.

As before, we assume that the data \(X(t)\) are observed with process
stochasticity around a mean function \(g(X(t))\)

\begin{equation}
X(t+1) = g(X(t)) + \varepsilon,
\end{equation}

where \(\varepsilon\) are IID normal random variables with zero-mean and
variance \(\sigma^2\). Note that we have chosen to assume additive
stochasticity. While we could consider log-normal stochasticity as in
the parametric models, we make this choice to emphasize that the
Gaussian process approach need not have structurally correct
stochasticity to be effective.

In order to make predictions, we update the GP based on the observed set
of transitions. To do so, we collect the time series of observed states
into a vector of ``current'' states,
\(\mathbf{X}_{\textrm{obs}} = \{X(1), \dots, X(T-1)\}\) and a vector of
``next'' states \(\mathbf{Y}_{\textrm{obs}} = \{X(2),\dots,X(T)\}\)
where \(T\) is the time of the final observation. Conditional on these
observations, the predicted next state, \(g(X_p)\), for any given
``current'' state, \(X_p\) follows a normal distribution with mean \(E\)
and variance \(C\) determined using the standard rules for conditioning
in multivariate normals, i.e.

\begin{equation}
E = K(X_p, \mathbf{X}_{\textrm{obs}}) \left(K(\mathbf{X}_{\textrm{obs}},\mathbf{X}_{\textrm{obs}}) - \sigma \mathbb{I}_n \right)^{-1} \mathbf{Y}_{\textrm{obs}}
\end{equation}

and

\begin{equation}
C = K(X_p, X_p) - K(X_p, \mathbf{X}_{\textrm{obs}}) \left(K(\mathbf{X}_{\textrm{obs}},\mathbf{X}_{\textrm{obs}}) - \sigma \mathbb{I} \right)^{-1} K(\mathbf{X}_{\textrm{obs}}, X_p)
\end{equation}

NA

\begin{equation}
K_{i,j}(\mathbf{x}, \mathbf{y}) = \exp\left( \frac{ -(x_i - y_j)^2}{2 \ell^2} \right)
\end{equation}

so that \(\ell\) gives the characteristic length-scale over which
correlation between two observations decays. See Rasmussen and Williams
(2006) for other choices of covariance kernels and their properties.
Note that this simple formulation assumes a prior mean of zero. For the
parameters we use inverse Gamma priors on both the length-scale \(\ell\)
and \(\sigma\), thus for example

\begin{equation}
f(\ell; \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} \ell^{-\alpha - 1}\exp\left(-\frac{\beta}{\ell}\right)
\end{equation}

For the prior on \(\ell\), \(\alpha =\) 10 and \(\beta =\) 10. The prior
on \(\sigma\), \(\alpha =\) 5 and \(\beta =\) 5.

We use a Metropolis-Hastings Markov Chain Monte Carlo (Gelman et al.
(2003)) to infer posterior distributions of the parameters of the GP
(Figure S4, code in appendix). Since the posterior distributions differ
substantially from the priors (Figure S4), most of the information in
the posterior comes from the data rather than the prior belief.

\subsection{The method of Gaussian Process Dynamic Programming
(GPDP)}\label{the-method-of-gaussian-process-dynamic-programming-gpdp}

We derive the harvest policy from the estimated GP by inserting it into
a SDP algorithm. Given the GP posteriors, we construct the transition
matrix representing the probability of going to each state \(X(t+1)\)
given any current state \(x(t)\) and any harvest \(h(t)\) (See the
function \texttt{gp\_transition\_matrix()} in the provided R package).
Given this transition matrix, we use the same value iteration algorithm
as in the parametric case to determine the optimal policy. In doing so,
the uncertainty in the future state under the GP, \(X(t+1)\), includes
both process uncertainty (based on the estimation of \(\sigma\)) and
structural uncertainty of the posterior collection of curves.

\section{Results}\label{results}

\subsection{Parametric and GP models for population
dynamics}\label{parametric-and-gp-models-for-population-dynamics}

To ensure our results are robust to the choice of parameters, we will
consider 96 different scenarios. To help better understand the process,
we first describe in detail the results of a single scenario.

\begin{figure}[htbp]
\centering
\includegraphics{components/figure/manuscript-figure_1-1.pdf}
\caption{Points show the training data of stock-size over time. Curves
show the expectations over the posterior step-ahead predictions based on
each of the estimated models. Observe that all models are fitting the
data reasonably well.)}
\end{figure}

All of the models fit the observed data rather closely and with
relatively small uncertainty. In Figure 1, we show the posterior
predictive curves. The training data of stock sizes observed over time
are points, overlaid with the step-ahead predictions of each estimated
model using the parameters sampled from their posterior distributions.
Compared to the true model most estimates appear to over-fit, predicting
patterns that are actually due purely to stochasticity. Model selection
criteria (Table 1) penalize more complex models and show a preference
for the simpler Ricker model over the models with alternative stable
states (Allen and Myers). Supplement provides details on the model
estimates.

\begin{longtable}[c]{@{}cccc@{}}
\caption{Model selection scores for several common criteria (DIC:
Deviance Information Criterion, AIC: Akaike Information Criterion, BIC:
Bayesian Information Criterion) all select the wrong 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.}\\
\toprule\addlinespace
\begin{minipage}[b]{0.12\columnwidth}\centering
~
\end{minipage} & \begin{minipage}[b]{0.10\columnwidth}\centering
Allen
\end{minipage} & \begin{minipage}[b]{0.11\columnwidth}\centering
Ricker
\end{minipage} & \begin{minipage}[b]{0.11\columnwidth}\centering
Myers
\end{minipage}
\\\addlinespace
\midrule\endhead
\begin{minipage}[t]{0.12\columnwidth}\centering
\textbf{DIC}
\end{minipage} & \begin{minipage}[t]{0.10\columnwidth}\centering
50.75
\end{minipage} & \begin{minipage}[t]{0.11\columnwidth}\centering
50.45
\end{minipage} & \begin{minipage}[t]{0.11\columnwidth}\centering
50.41
\end{minipage}
\\\addlinespace
\begin{minipage}[t]{0.12\columnwidth}\centering
\textbf{AIC}
\end{minipage} & \begin{minipage}[t]{0.10\columnwidth}\centering
-24.51
\end{minipage} & \begin{minipage}[t]{0.11\columnwidth}\centering
-30.13
\end{minipage} & \begin{minipage}[t]{0.11\columnwidth}\centering
-27.01
\end{minipage}
\\\addlinespace
\begin{minipage}[t]{0.12\columnwidth}\centering
\textbf{BIC}
\end{minipage} & \begin{minipage}[t]{0.10\columnwidth}\centering
-17.75
\end{minipage} & \begin{minipage}[t]{0.11\columnwidth}\centering
-25.06
\end{minipage} & \begin{minipage}[t]{0.11\columnwidth}\centering
-20.25
\end{minipage}
\\\addlinespace
\bottomrule
\end{longtable}

\begin{figure}[htbp]
\centering
\includegraphics{manuscript-figure_2-1.jpeg}
\caption{The inferred Gaussian process compared to the true process and
maximum-likelihood estimated (MLE) process. We show the expected value
for the function \(f\) under each model. Plots show replicates drawn
from the posterior distributions in order to convey uncertainty of the
estimates. Note the MLE is a point estimate of parameters and so
reflects no uncertainty in the distribution. The training data are also
shown as black points. The GP is conditioned on (0,0), shown as a
pseudo-data point.}
\end{figure}

We show the mean inferred population dynamics of each model relative to
the true model used to generate the data in Figure 2, predicting the
relationship between observed population size (x-axis) to the population
size after recruitment the following year.\\In addition to the raw data,
the GP is conditioned on going through the point 0,0 without error. All
parametric models also make this assumption. Conditioning on (0,0) is
equivalent to making the assumption that the population is closed, so
that once it hits 0 it stays at 0, despite the lack of any data in the
observed sequence to justify this. This assumption illustrates how the
GP can capture common-sense biology without having to assume more
explicit details about the dynamics at low population numbers that have
never been observed. If the population were not closed, one could repeat
the entire analysis without this assumption. Unlike parametric models,
the GP corresponds to a distribution of curves, of which this plot only
shows the means. Uncertainty in the parameters of the GP (not shown)
further widens the band of possible population sizes. In Figure S1 (see
supplement), we show the performance of the models outside the observed
training data.

\begin{figure}[htbp]
\centering
\includegraphics{components/figure/manuscript-figure_3-1.pdf}
\caption{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 (Reed 1979) Several policies show a numerical
jitter due to the discretization of states in the dynamic programming
algorithm; doubling the number of grid points did not qualitatively
change the results.}
\end{figure}

Despite the similarities in model fits to the observed data, the
policies inferred under each model differ widely (Figure 3). Policies
are shown in terms of target escapement, \(S(t) = x_t - h\). Under
models such as this a constant escapement policy is expected to be
optimal (Reed 1979), whereby population levels below a certain size
\(S\) are unharvested, while above that size the harvest strategy aims
to return the population to \(S\). Whenever a model predicts that the
population will not persist below a certain threshold, the optimal
solution is to harvest the entire population immediately, resulting in
an escapement \(S=0\), as seen in the true (correct form, exact
parameters) model, the Allen model (correct form, estimated parameters)
and the GP. Only the structurally correct model (Allen model) and the GP
produce policies close to the true optimum policy.

\begin{figure}[htbp]
\centering
\includegraphics{manuscript-figure_4-1.jpeg}
\caption{In the management context, GPDP outperforms approaches based on
parametric models. We show 100 replicate simulations of the stock
dynamics (Eqn 1) under the policies derived from each of the estimated
models, as well as the policy based on the exact underlying model.}
\end{figure}

In Figure 4, we show the consequences of managing 100 replicate
realizations of the simulated fishery under policies derived from each
model. The structurally correct model under-harvests, leaving the stock
to vary around its unfished optimum. The structurally incorrect Ricker
model over-harvests the population past the tipping point consistently,
resulting in the immediate crash of the stock and thus leads to minimal
long-term catch.

The results across replicate stochastic simulations are most easily
compared by using the relative differences in net present value realized
by each of the model (Figure 5). Although not perfect, the GPDP
consistently realizes a value close to the optimal solution, and avoids
ever driving the system across the tipping point, which results in the
near-zero value cases in the parametric models.

\begin{figure}[htbp]
\centering
\includegraphics{components/figure/manuscript-figure_5-1.pdf}
\caption{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.}
\end{figure}

\subsection{Sensitivity Analysis}\label{sensitivity-analysis}

These results are not sensitive to the modeling details of the
simulation. The GPDP estimate remains very close to the optimal solution
(obtained by knowing the true model) across changes to the training
simulation, scale of stochasticity, parameters or structure of the
underlying model. In the Supplement, we consider both a Latin hypercube
approach and a more focused investigation of the effects of the relative
distance to the Allee threshold and the variance of process
stochasticity.

The GPDP is only weakly influenced by increasing stochasticity or
increasing Allee effects over much of the range (Figure S2). Larger
\(\sigma\) or higher Allee levels make even the optimal solution without
any model or parameter uncertainty unable to harvest the population
effectively (e.g.~the stochasticity is large enough to violate the
self-sustaining criterion of Reed (1979)).

\section{Discussion}\label{discussion}

Simple, mechanistically motivated models offer the potential to increase
our basic understanding of ecological processes (Geritz and Kisdi 2012,
Cuddington et al. 2013). But such models can be both inaccurate and
misleading when used in a decision making framework. In this paper we
tackled two aspects of uncertainty that are common to many ecological
decision-making problems and fundamentally challenging to existing
approaches that largely rely on parametric models:

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\itemsep1pt\parskip0pt\parsep0pt
\item
  We do not know the correct models for ecological systems.
\item
  We have limited data from which to estimate the model.
\end{enumerate}

We have illustrated how the use of non-parametric methods provides more
reliable solutions in the sequential decision-making problem.

\subsubsection{Traditional model-choice approaches can be positively
misleading.}\label{traditional-model-choice-approaches-can-be-positively-misleading.}

Our results illustrate that model-choice approaches can be absolutely NA
dynamics because they have fewer parameters and the data are far from
the tipping point. That is, when the data come from around the stable
steady state, all the parametric models are approximately linear and
approximately identical. Thus, it is intuitive that all model selection
methods choose the simplest model. In a complex world, the result is
suboptimal. But in a world that might contain tipping points, the result
could be disastrous.

Many managers both in fisheries and beyond face a similar situation:
they have not observed the population dynamics at all possible
densities. A lack of comprehensive data at all population sizes,
combined with the inability to formulate accurate population models for
low population sizes in the absence of data, makes this situation the
rule more than the exception. Relying on parametric models and model
choice processes that favor simplicity ignores this basic reality. For a
long time, Carl Walters (e.g. Walters and Hilborn 1978) has argued that
if we began by fishing any newly exploited population down to very low
levels and then let it recover, we would be much better at estimating
population dynamics and thus predicting the optimal harvest levels.
While certainly true, this presents a rather risky policy in the face of
potential tipping points. The GPDP offers a risk-adverse alternative.

\subsubsection{GPDP population dynamics capture larger uncertainty in
regions where the data are
poor}\label{gpdp-population-dynamics-capture-larger-uncertainty-in-regions-where-the-data-are-poor}

Parametric models perform most poorly when we seek a management strategy
outside the range of the observed data. The GPDP, in contrast, leads to
a predictive model that expresses a great deal of uncertainty about the
probable dynamics \emph{outside} the range of the observed data, while
retaining very good predictive accuracy \emph{inside} the range. The
management policy based on by the GPDP balances uncertainty outside the
range of the observed data against the immediate value of the harvest,
and acts to stabilize the population dynamics in a region of state space
in which the predictions are reliably reflected by the data.

Such problems are ubiquitous across ecological decision-making and
conservation where the greatest concerns involve decisions that lead to
population sizes that have never been observed and for which we do not
know the response -- whether this is the collapse of a fishery, the
spread of an invasive, or the loss of habitat.

\subsubsection{The role of the prior}\label{the-role-of-the-prior}

Outside of the observed range of the data, the GP reverts to the prior,
and consequently the choice of the prior can also play a significant
role in determining the optimal policy. In the examples shown here we
have selected a prior that is both relatively uninformative (due to the
broad priors placed on its parameters \(\ell\) and \(\sigma\) and simple
(the choice of our covariance function, Eqns 12 and 13 ). In practice,
these should be chosen to confer particular biological properties. In
principle, this may allow a manager to improve the performance of the
GPDP by adding detail as is justified. For instance, it would be
possible to use a linear or a Ricker-shaped mean in the prior without
making the much stronger assumption that the Ricker is the structurally
correct model (Sugeno and Munch 2013a). One fruitful avenue of future
research is identifying criteria to ensure the prior and the reward
function are chosen appropriately for the problem at hand.

\section{Acknowledgments}\label{acknowledgments}

This work was partially supported by NOAA-IAM grant to SM and Alec
McCall and administered through the Center for Stock Assessment
Research, a partnership between the University of California Santa Cruz
and the Fisheries Ecology Division, Southwest Fisheries Science Center,
Santa Cruz, CA and by NSF grant EF-0924195 to MM and NSF grant
DBI-1306697 to CB.

Allen, D., and K. Tanner. 2005. Infusing active learning into the
large-enrollment biology class: seven strategies, from the simple to
complex. Cell biology education 4:262--8.

Athanassoglou, S., and A. Xepapadeas. 2012. Pollution control with
uncertain stock dynamics: When, and how, to be precautious. Journal of
Environmental Economics and Management 63:304--320.

Boettiger, C., M. Mangel, and S. Munch. 2014. nonparametric-bayes
v0.1.0. ZENODO. \url{http://dx.doi.org/10.5281/zenodo.12669}.

Brozović, N., and W. Schlenker. 2011. Optimal management of an ecosystem
with an unknown threshold. Ecological Economics:1--14.

Burnham, K. P., and D. R. Anderson. 2002. Model Selection and
Multi-Model Inference. Page 496. Springer.

Clark, C. W. 1976. Mathematical Bioeconomics. WileyNew York.

Clark, C. W. 2009. Mathematical Bioeconomics. WileyNew York.

Clark, C. W., and G. P. Kirkwood. 1986. On uncertain renewable resource
stocks: Optimal harvest policies and the value of stock surveys. Journal
of Environmental Economics and Management 13:235--244.

Clark, C. W., and M. Mangel. 2000. Dynamic state variable models in
ecology. Oxford University PressOxford.

Courchamp, F., L. Berec, and J. Gascoigne. 2008. Allee Effects in
Ecology and Conservation. Page 256. Oxford University Press, USA.

Cressie, N., C. a Calder, J. S. Clark, J. M. {Ver Hoef}, and C. K.
Wikle. 2009. Accounting for uncertainty in ecological analysis: the
strengths and limitations of hierarchical statistical modeling.
Ecological Applications 19:553--70.

Cuddington, K. M., M. Fortin, and L. Gerber. 2013. Process-based models
are required to manage ecological systems in a changing world. Ecosphere
4:1--12.

Fischer, J., G. D. Peterson, T. a Gardner, L. J. Gordon, I. Fazey, T.
Elmqvist, A. Felton, C. Folke, and S. Dovers. 2009. Integrating
resilience thinking and optimisation for conservation. Trends in ecology
\& evolution 24:549--54.

Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin. 2003. Bayesian
Data Analysis. 2nd editions. Chapman; Hall/CRC.

Geritz, S. A. H., and E. Kisdi. 2012. Mathematical ecology: why
mechanistic models? Journal of mathematical biology 65:1411--5.

Gordon, H. 1954. The economic theory of a common-property resource: the
fishery. The Journal of Political Economy 62:124--142.

Hilborn, R., and M. Mangel. 1997. The Ecological Detective: Confronting
Models with data. Page 330. Princeton University Press.

Hughes, T. P., C. Linares, V. Dakos, I. a van de Leemput, and E. H. van
Nes. 2013. Living dangerously on borrowed time during slow, unrecognized
regime shifts. Trends in ecology \& evolution 28:149--55.

Kocijan, J., A. Girard, B. Banko, and R. Murray-Smith. 2005. Dynamic
systems identification with Gaussian processes. Mathematical and
Computer Modelling of Dynamical Systems 11:411--424.

Ludwig, D., and C. J. Walters. 1982. Optimal harvesting with imprecise
parameter estimates. Ecological Modelling 14:273--292.

Mangel, M. 2014. Stochastic Dynamic Programming Illuminates the Link
Between Environment. Bulletin of Mathematical Biology in press.

Mangel, M., and C. W. Clark. 1988. Dynamic Modeling in Behavioral
Ecology. (J. Krebs and T. Clutton-Brock, Eds.). Princeton University
PressPrinceton.

Mangel, M., 0. Fiksen, and J. Giske. 2001. Theoretical and statistical
models in natural resource management and research. Pages 57--71
\emph{in} T. M. Shenk and A. B. Franklin, editors. Modeling in natural
resource management, development, interpretation and application. Island
PressWashington DC.

Marescot, L., G. Chapron, I. Chadès, P. L. Fackler, C. Duchamp, E.
Marboutin, and O. Gimenez. 2013. Complex decisions made simple: a primer
on stochastic dynamic programming. Methods in Ecology and
Evolution:n/a--n/a.

May, R. M., J. R. Beddington, C. W. Clark, S. J. Holt, and R. M. Laws.
1979. Management of multispecies fisheries. Science (New York, N.Y.)
205:267--77.

McAllister, M. 1998. Bayesian stock assessment: a review and example
application using the logistic model. ICES Journal of Marine Science
55:1031--1060.

Munch, S. B., A. Kottas, and M. Mangel. 2005a. Bayesian nonparametric
analysis of stock-recruitment relationships. Canadian Journal of
Fisheries and Aquatic Sciences 62:1808--1821.

Munch, S. B., M. L. Snover, G. M. Watters, and M. Mangel. 2005b. A
unified treatment of top-down and bottom-up control of reproduction in
populations. Ecology Letters 8:691--695.

Myers, R. A., N. J. Barrowman, J. A. Hutchings, and A. a Rosenberg.
1995. Population dynamics of exploited fish stocks at low population
levels. Science (New York, N.Y.) 269:1106--8.

Polasky, S., S. R. Carpenter, C. Folke, and B. Keeler. 2011.
Decision-making under great uncertainty: environmental management in an
era of global change. Trends in ecology \& evolution:1--7.

R Core Team. 2013. R: A Language and Environment for Statistical
Computing. R Foundation for Statistical ComputingVienna, Austria.

Rasmussen, C. E., and C. K. I. Williams. 2006. Gaussian Processes for
Machine Learning. (Thomas Dietterich, Ed.). MIT Press,Boston.

Reed, W. J. 1979. Optimal escapement levels in stochastic and
deterministic harvesting models. Journal of Environmental Economics and
Management 6:350--363.

Roughgarden, J. E., and F. Smith. 1996. Why fisheries collapse and what
to do about it. Proceedings of the National Academy of Sciences of the
United States of America 93:5078.

Schapaugh, A. W., and A. J. Tyre. 2013. Accounting for parametric
uncertainty in Markov decision processes. Ecological Modelling
254:15--21.

Scheffer, M., J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V.
Dakos, H. Held, E. H. van Nes, M. Rietkerk, and G. Sugihara. 2009.
Early-warning signals for critical transitions. Nature 461:53--9.

Scheffer, M., S. R. Carpenter, J. A. Foley, C. Folke, and B. Walker.
2001. Catastrophic shifts in ecosystems. Nature 413:591--6.

Sethi, G., C. Costello, A. Fisher, M. Hanemann, and L. Karp. 2005.
Fishery management under multiple uncertainty. Journal of Environmental
Economics and Management 50:300--318.

Sigourney, D. B., S. B. Munch, and B. H. Letcher. 2012. Combining a
Bayesian nonparametric method with a hierarchical framework to estimate
individual and temporal variation in growth. Ecological Modelling
247:125--134.

Su, Y.-S., and Masanao Yajima. 2013. R2jags: A Package for Running jags
from R.

Sugeno, M., and S. B. Munch. 2013a. A semiparametric Bayesian method for
detecting Allee effects. Ecology 94:1196--1204.

Sugeno, M., and S. B. Munch. 2013b. A semiparametric Bayesian approach
to estimating maximum reproductive rates at low population sizes.
Ecological applications : a publication of the Ecological Society of
America 23:699--709.

Thorson, J. T., K. Ono, and S. B. Munch. 2014. A Bayesian approach to
identifying and compensating for model misspecification in population
models. Ecology 95:329--41.

Walters, C. J., and R. Hilborn. 1978. Ecological Optimization and
Adaptive Management. Annual Review of Ecology and Systematics
9:157--188.

Weitzman, M. L. 2013. A Precautionary Tale of Uncertain Tail Fattening.
Environmental and Resource Economics 55:159--173.

Williams, B. K. 2001. Uncertainty , learning , and the optimal
management of wildlife. Environmental and Ecological Statistics
8:269--288.

Xie, Y. 2013. Dynamic documents with R and knitr. Chapman; Hall/CRCBoca
Raton, Florida.

\end{document}

repo_name: cboettig/nonparametric-bayes

repo_language: TeX

mime_type: text/plain


text:
FROM rocker/r-ver:3.2.0
MAINTAINER Carl Boettiger cboettig@ropensci.org

## Install software dependencies
RUN apt-get update \
  && apt-get install -y \
    libssl-dev \
    libxml2-dev \
    libcurl4-openssl-dev \
    libnlopt-dev \
    libpng-dev \
    git \
    wget \
#    pandoc \
#    pandoc-citeproc \
    libgsl0-dev \
    libgl1-mesa-dev \
    libglu1-mesa-dev
# wget https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/Source/JAGS-4.2.0.tar.gz -O jags.tar.gz \
RUN wget https://sourceforge.net/projects/mcmc-jags/files/JAGS/3.x/Source/JAGS-3.4.0.tar.gz/download  -O jags.tar.gz \
  && tar -xf jags.tar.gz \
  && cd JAGS* && ./configure && make && make install \
  && cd / && rm -rf jags.tar.gz JAGS*
RUN . /etc/environment \
  && R -e 'install.packages(c("devtools")); devtools::install_github("cboettig/pdg_control")' \
  && R -e 'devtools::install_github(c("cboettig/nonparametric-bayes", "cboettig/cboettigR@8bb11be44195e70dfd731cd98f11db1b73e62303"), dep=TRUE)'

## Not needed for code, but only to compile Rmd to pdf:
RUN apt-get update && apt-get -y install \
  texlive-latex-recommended texlive-fonts-recommended

## debian:jessie version of pandoc-citeproc is too old.
## also note template is not compatible with newer versions of pandoc.
RUN apt-get update && apt-get -y install \
  cabal-install libghc-cabal-dev cpphs libghc-pretty-show-dev libexpat1-dev
RUN cabal update && cabal install hsb2hs
RUN cabal install pandoc-1.13 --force-reinstalls
RUN cabal install pandoc-citeproc-0.5
RUN ln -s /root/.cabal/bin/pandoc /usr/bin/pandoc
RUN ln -s /root/.cabal/bin/pandoc-citeproc /usr/bin/pandoc-citeproc

COPY . /nonparametric-bayes
WORKDIR /nonparametric-bayes


text:
Sources for: Avoiding tipping points in fisheries management through Gaussian Process Dynamic Programming
========================================================================================================

[![Build Status](http://server.carlboettiger.info:88/api/badge/github.com/cboettig/nonparametric-bayes/status.svg?branch=master)](http://server.carlboettiger.info:88/github.com/cboettig/nonparametric-bayes)

<!--
Manuscript code is run automatically any time changes are made to
this repository.  The badge above summarizes the status and contains
links to further details.
-->

This directory contains the `.Rmd` files for the manuscript and
supplement. These files use functions from the R package from this
repository (`nonparametric-bayes`) but also several additional R packages
and other software which are responsible for the parametric comparisons,
as well as software responsible merely for the formatting to go from
`.Rmd` to a `.pdf` file (e.g. `knitr`, `pandoc`, `LaTeX`).

All of these additional dependencies means that rebuilding the
manuscript from scratch means more software to install and more
points of failure. To help you get started, we provide a copy
of the software environment we use in this paper as a Docker
container. If you'd rather not worry about the installation details,
[Install Docker](https://docs.docker.com/installation)
on your laptop or a [user-friendly cloud
server](https://www.digitalocean.com/?refcode=08c6ac401e49) [^1]

The badge above is automatically generated when the manuscript
is rebuilt on a remote server. This provides a convenient way
to confirm that everything is still working.

Docker Quickstart
-----------------

[Install Docker](https://docs.docker.com/installation) on your laptop or server. Then from the `manuscripts/` directory run:

```bash
docker build -t cboettig/nonparametric-bayes .
```

```
docker run --rm -ti -v $(pwd):/data cboettig/nonparametric-bayes R -e "rmarkdown::render('manuscript.Rmd')"
```

You can now open `manuscript.Rmd` or `supplement.Rmd` files in RStudio and run the code interactively or compile the pdfs from scratch.  See [RStudio's rmarkdown](http://rmarkdown.rstudio.com/) for details.

Local Installation requirements
--------------------------------

- RStudio (includes `pandoc`)
- A LaTeX environment
- The `nonparametric-bayes` package [from github](http://github.com/cboettig/nonparametric-bayes), including all packages in the [DESCRIPTION](../DESCRIPTION) file (by default suggested packages are not installed).

--------

[1]: this link provides $10 credit to Digital Ocean; enough to run a 1GB instance for a month.  Disclosure: If you subsequently become a Digital Ocean customer the authors receive a credit for the referral.


text:
#!/usr/bin/R

# restore cache file
archive <- "nonparametric-bayes-cache.tar.gz"

setwd("components")
download.file(paste0("http://two.ucdavis.edu/~cboettig/data/", archive), archive)
untar(archive)
unlink(archive)
setwd("..")


text:
sensitivity <- function(model, parameters, nuisance, seed){

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

sigma_g <- nuisance[["sigma_g"]]
p <- parameters
sigma_m <- 0.0
z_g <- function() rlnorm(1, 0, sigma_g)
z_m <- function() 1
x_grid <- seq(0, 15, length=50)
h_grid <- x_grid
profit <- function(x,h) pmin(x, h)
delta <- 0.01
OptTime <- 50  # stationarity with unstable models is tricky thing
reward <- 0
xT <- 0
Xo <- 5.5# observations start from
x0 <- 8 # simulation under policy starts from
Tobs <- 40
MaxT <- 1000 # timeout for value iteration convergence

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

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

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

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

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

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

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

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

                     sims <- lapply(OPT, function(D){
                       set.seed(1)
                       lapply(1:100, function(i)
                         ForwardSimulate(f, p, x_grid, h_grid, x0, D, z_g, profit=profit, OptTime=OptTime)
                              )
                       })

                     dat <- melt(sims, id=names(sims[[1]][[1]]))
                     sims_data <- data.table(dat)
                     setnames(sims_data, c("L1", "L2"), c("method", "reps"))
                     # Legend in original ordering please, not alphabetical:
                     sims_data$method = factor(sims_data$method, ordered=TRUE, levels=names(OPT))

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

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

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

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

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


text:
<?xml version="1.0" encoding="utf-8"?>
<style xmlns="http://purl.org/net/xbiblio/csl" class="in-text" default-locale="en-US" version="1.0" demote-non-dropping-particle="sort-only">
  <info>
    <title>Ecology</title>
    <id>http://www.zotero.org/styles/ecology</id>
    <link href="http://www.zotero.org/styles/ecology" rel="self"/>
    <link href="http://esapubs.org/esapubs/AuthorInstructions.htm" rel="documentation"/>
    <author>
      <name>Rintze Zelle</name>
      <uri>http://twitter.com/rintzezelle</uri>
    </author>
    <category citation-format="author-date"/>
    <category field="biology"/>
    <issn>0012-9658</issn>
    <updated>2012-09-27T22:06:38+00:00</updated>
    <rights license="http://creativecommons.org/licenses/by-sa/3.0/">This work is licensed under a Creative Commons Attribution-ShareAlike 3.0 License</rights>
  </info>
  <macro name="container-contributors">
    <choose>
      <if type="chapter paper-conference" match="any">
        <text term="in" suffix=" " font-style="italic"/>
        <names variable="editor translator" delimiter=", ">
          <name and="text" initialize-with=". " delimiter=", "/>
          <label form="long" prefix=", " suffix="."/>
        </names>
      </if>
    </choose>
  </macro>
  <macro name="secondary-contributors">
    <choose>
      <if type="chapter paper-conference" match="none">
        <names variable="editor translator" delimiter=", " prefix=" (" suffix=")">
          <name and="text" initialize-with=". " delimiter=", "/>
          <label form="short" prefix=", " text-case="capitalize-first"/>
        </names>
      </if>
    </choose>
  </macro>
  <macro name="author">
    <names variable="author">
      <name name-as-sort-order="first" and="text" sort-separator=", " initialize-with=". " delimiter=", " delimiter-precedes-last="always"/>
      <label form="short" prefix=" (" suffix=".)" text-case="capitalize-first" strip-periods="true"/>
      <substitute>
        <names variable="editor"/>
        <names variable="translator"/>
        <text variable="title"/>
      </substitute>
    </names>
  </macro>
  <macro name="author-short">
    <names variable="author">
      <name form="short" and="text" delimiter=", " initialize-with=". "/>
      <substitute>
        <names variable="editor"/>
        <names variable="translator"/>
        <choose>
          <if type="bill book graphic legal_case legislation motion_picture report song" match="any">
            <text variable="title" form="short" font-style="italic"/>
          </if>
          <else>
            <text variable="title" form="short" quotes="true"/>
          </else>
        </choose>
      </substitute>
    </names>
  </macro>
  <macro name="access">
    <choose>
      <if type="webpage">
        <text variable="URL"/>
      </if>
    </choose>
  </macro>
  <macro name="publisher">
    <group delimiter=", ">
      <choose>
        <if type="article-journal article-magazine" match="none">
          <text variable="genre"/>
          <text variable="publisher"/>
          <text variable="publisher-place"/>
        </if>
      </choose>
    </group>
  </macro>
  <macro name="event">
    <text variable="event"/>
    <text variable="event-place" prefix=", "/>
  </macro>
  <macro name="issued">
    <choose>
      <if variable="issued">
        <group prefix=" " suffix=".">
          <date variable="issued">
            <date-part name="year"/>
          </date>
          <choose>
            <if type="article-journal bill book chapter graphic legal_case legislation motion_picture paper-conference report song" match="none">
              <date variable="issued">
                <date-part prefix=", " name="month"/>
                <date-part prefix=" " name="day"/>
              </date>
            </if>
          </choose>
        </group>
      </if>
      <else>
        <text prefix=" (" term="no date" suffix=")." form="short"/>
      </else>
    </choose>
  </macro>
  <macro name="issued-year">
    <choose>
      <if variable="issued">
        <date variable="issued">
          <date-part name="year"/>
        </date>
      </if>
      <else>
        <text term="no date" form="short"/>
      </else>
    </choose>
  </macro>
  <macro name="edition">
    <choose>
      <if type="bill book chapter graphic legal_case legislation motion_picture paper-conference report song" match="any">
        <choose>
          <if is-numeric="edition">
            <number variable="edition" form="long-ordinal" text-case="capitalize-first"/>
            <text term="edition" form="long" prefix=" " suffix="."/>
          </if>
          <else>
            <text variable="edition" suffix="."/>
          </else>
        </choose>
      </if>
    </choose>
  </macro>
  <macro name="locators">
    <choose>
      <if type="article-journal article-magazine article-newspaper" match="any">
        <text variable="container-title" prefix=". "/>
        <text variable="volume" prefix=" "/>
        <text variable="page" prefix=":"/>
      </if>
      <else-if type="bill book chapter graphic legal_case legislation motion_picture paper-conference report song" match="any">
        <group prefix=". " delimiter=" ">
          <label variable="page" form="long" text-case="capitalize-first"/>
          <text variable="page"/>
          <text macro="container-contributors"/>
          <text macro="secondary-contributors"/>
          <text variable="container-title"/>
        </group>
      </else-if>
    </choose>
  </macro>
  <macro name="citation-locator">
    <group>
      <label variable="locator" form="short"/>
      <text variable="locator" prefix=" "/>
    </group>
  </macro>
  <citation et-al-min="3" et-al-use-first="1" disambiguate-add-year-suffix="true" collapse="year">
    <sort>
      <key variable="issued"/>
      <key macro="author"/>
    </sort>
    <layout prefix="(" suffix=")" delimiter=", ">
      <group delimiter=" ">
        <text macro="author-short"/>
        <text macro="issued-year"/>
        <text macro="citation-locator"/>
      </group>
    </layout>
  </citation>
  <bibliography hanging-indent="true" entry-spacing="0">
    <sort>
      <key macro="author"/>
      <key variable="issued"/>
    </sort>
    <layout suffix=".">
      <text macro="author" suffix="."/>
      <text macro="issued" suffix=" "/>
      <text variable="title"/>
      <text macro="locators"/>
      <group delimiter=". " prefix=". ">
        <text macro="edition"/>
        <text macro="publisher"/>
        <text macro="access"/>
      </group>
    </layout>
  </bibliography>
</style>


text:
t: latex
bibliography: components/references.bib
csl: components/ecology.csl
template: components/elsarticle.latex
latex-engine: pdflatex


text:
## @knitr plotting-options
library(knitr)
source("~/.knitr_defaults.R")
#opts_knit$set(upload.fun = socialR::flickr.url)
library(knitcitations)
library(nonparametricbayes)

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

## @knitr stateeq
f <- RickerAllee
p <- c(2, 8, 5)
K <- 10  # approx, a li'l' less
allee <- 5 # approx, a li'l' less

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

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

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

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

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

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

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

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

## @knitr gp-output
# Summarize the GP model
tgp_dat <-
    data.frame(  x = x_grid,
                 y = gp_dat$E_Ef,
                 ymin = gp_dat$E_Ef - 2 * sqrt(gp_dat$E_Vf),
                 ymax = gp_dat$E_Ef + 2 * sqrt(gp_dat$E_Vf) )

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

## @knitr allen-posteriors
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)

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

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

## @knitr ricker-posteriors
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)

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

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

## @knitr myers-posteriors
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)

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

## @knitr assemble-models
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
require(MASS)
step_ahead <- function(x, f, p){
  h = 0
  x_predict <- sapply(x, f, h, p)
  n <- length(x_predict) - 1
  y <- c(x[1], x_predict[1:n])
  y
}

## @knitr Figureb
names(bayes_pars) <- c("r0", "K", "theta")
names(myers_bayes_pars) <- c("r0",  "theta", "K")
names(ricker_bayes_pars) <- c("r0",  "K")
gp_f_at_obs <- gp_predict(gp, x, burnin=1e4, thin=300)
df <- melt(data.frame(time = 1:length(x), stock = x,
                GP = gp_f_at_obs$E_Ef,
                True = step_ahead(x,f,p),
                MLE = step_ahead(x,f,est$p),
                Allen = step_ahead(x, allen_f, bayes_pars),
                Ricker = step_ahead(x, ricker_f, ricker_bayes_pars),
                Myers = step_ahead(x, myers_f, myers_bayes_pars)
                 ), id=c("time", "stock"))
Figure1b <- ggplot(df) + geom_point(aes(time, stock)) +
  geom_line(aes(time, value, col=variable)) +
    scale_colour_manual(values=colorkey)

## @knitr Figureb_posteriors
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

## @knitr statespace_plot
statespace_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=1, alpha = 0.7) +
    geom_point(data = obs, aes(x,y), alpha = 0.8) +
    xlab(expression(X[t])) + ylab(expression(X[t+1])) +
    scale_colour_manual(values = colorkey)
statespace_plot

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

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

## @knitr gp-opt
# uses expected values from GP, instead of integrating over posterior
#matrices_gp <- gp_transition_matrix(gp_dat$E_Ef, gp_dat$E_Vf, x_grid, h_grid)
matrices_gp <- gp_transition_matrix(gp_dat$Ef_posterior, gp_dat$Vf_posterior, x_grid, h_grid)
opt_gp <- value_iteration(matrices_gp, x_grid, h_grid, MaxT, xT, profit, delta, reward)

## @knitr mle-opt
matrices_true <- f_transition_matrix(f, p, x_grid, h_grid, sigma_g)
opt_true <- value_iteration(matrices_true, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
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)

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

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

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

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

## @knitr Figure2
policies <- melt(data.frame(stock=x_grid, sapply(OPT, function(x) x_grid[x])), id="stock")
names(policies) <- c("stock", "method", "value")

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

## @knitr sims
sims <- lapply(OPT, function(D){
  set.seed(1)
  lapply(1:100, function(i)
    ForwardSimulate(f, p, x_grid, h_grid, x0, D, z_g, profit=profit, OptTime=OptTime)
  )
})

dat <- melt(sims, id=names(sims[[1]][[1]]))
sims_data <- data.table(dat)
setnames(sims_data, c("L1", "L2"), c("method", "reps"))
# Legend in original ordering please, not alphabetical:
sims_data$method = factor(sims_data$method, ordered=TRUE, levels=names(OPT))

## @knitr Figure3
ggplot(sims_data) +
  geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=.1) +
  scale_colour_manual(values=colorkey, guide = guide_legend(override.aes = list(alpha = 1)))

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

## @knitr Figure4plots
fig4v1 <- ggplot(actual_over_optimal, aes(value)) + geom_histogram(aes(fill=variable)) +
  facet_wrap(~variable, scales = "free_y")  + guides(legend.position = "none") +
  xlab("Total profit by replicate") + scale_fill_manual(values=colorkey) # density plots fail when delta fn
fig4v2 <- ggplot(actual_over_optimal, aes(value)) + geom_histogram(aes(fill=variable), binwidth=0.1) +
  xlab("Total profit by replicate")+ scale_fill_manual(values=colorkey)
fig4v3 <- 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)
fig4v3

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

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

## @knitr appendixplots
gp_assessment_plots[[1]]
gp_assessment_plots[[2]]
plot_allen_traces
plot_allen_posteriors
plot_ricker_traces
plot_ricker_posteriors
plot_myers_traces
plot_myers_posteriors