In cboettig/nonparametric-bayes: Nonparametric Bayes inference for ecological models

Comparison of Nonparametric Bayesian Gaussian Process estimates to standard the Parametric Bayesian approach

Plotting and knitr options, (can generally be ignored)

source("~/.knitr_defaults.R")
#opts_knit$set(upload.fun = socialR::flickr.url)
library(knitcitations)
library(nonparametricbayes) 
opts_chunk$set(external=TRUE)
read_chunk("external-chunks.R")

Model and parameters

f <- Myers
p <- c(1, 2, 6)
K <- 5  # approx, a li'l' less
allee <- 1.2 # approx, a li'l' less

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

Sample Data

Maximum Likelihood

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,.2), 
         p[2]*rlnorm(1,0,.1), 
         p[3]*rlnorm(1,0, .1), 
         sigma_g * rlnorm(1,0,.3))
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)

Mean predictions

Non-parametric Bayes

Estimate the Gaussian Process (nonparametric Bayesian fit)

Show traces and posteriors against priors

Parametric Bayesian Models

We use the JAGS Gibbs sampler, a recent open source BUGS implementation with an R interface that works on most platforms. We initialize the usual MCMC parameters; see ?jags for details.

All parametric Bayesian estimates use the following basic parameters for the JAGS MCMC:

We will use the same priors for process and observation noise in each model,

Parametric Bayes of correct (Allen) model

We initiate the MCMC chain (init_p) using the true values of the parameters p from the simulation. While impossible in real data, this gives the parametric Bayesian approach the best chance at succeeding. y is the timeseries (recall obs has the $x_t$, $x_{t+1}$ pairs)

The actual model is defined in a model.file that contains an R function that is automatically translated into BUGS code by R2WinBUGS. The file defines the priors and the model. We write the file from R as follows:

Write the priors into a list for later reference

We define which parameters to keep track of, and set the initial values of parameters in the transformed space used by the MCMC. We use logarithms to maintain strictly positive values of parameters where appropriate.

Convergence diagnostics for Allen model

R notes: this strips classes from the mcmc.list object (so that we have list of matrices; objects that reshape2::melt can handle intelligently), and then combines chains into one array. In this array each parameter is given its value at each sample from the posterior (index) for each chain.

Reshape the posterior parameter distribution data, transform back into original space, and calculate the mean parameters and mean function

Parametric Bayes based on the structurally wrong model (Ricker)

Compute prior curves

We define which parameters to keep track of, and set the initial values of parameters in the transformed space used by the MCMC.

Convergence diagnostics for parametric bayes Ricker model

Reshape posteriors data, transform back, calculate mode and corresponding function.

Myers Parametric Bayes

jags.params=c("r0", "theta", "K", "stdQ")
jags.inits <- function(){
  list("r0"= rlnorm(1,0,.1), 
       "K"=    6 * rlnorm(1,0,.1),
       "theta" = 2 * rlnorm(1,0,.1),  
       "stdQ"= 0.1 * rlnorm(1,0,.1),
       .RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
set.seed(12345)
myers_jags <- do.call(jags.parallel, 
                      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"))

Convergence diagnostics for parametric bayes

Phase-space diagram of the expected dynamics

Goodness of fit

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

Optimal policies by value iteration

Compute the optimal policy under each model using stochastic dynamic programming. We begin with the policy based on the GP model,

Determine the optimal policy based on the allen and MLE models

Determine the optimal policy based on Bayesian Allen model

Bayesian Ricker

Bayesian Myers model

Assemble the data

Graph of the optimal policies

Simulate 100 realizations managed under each of the policies

cboettig/nonparametric-bayes documentation built on May 13, 2019, 2:09 p.m.