Computation time is an unavoidable consideration when working with simulation-based inference, for all but small datasets and simple models.
The pomp package therefore allows you to specify the most computationally intensive steps—usually, simulation of the stochastic dynamic system, and evaluation of the measurement density—as snippets of C code.
Consequently, to use pomp, your R program must have access to a C compiler. In addition, pomp takes advantage of some Fortran code and therefore requires a Fortran compiler.
Installing the necessary compilers should be fairly routine, but does involve an extra step beyond the usual installation of an R package, unless you are running the Linux operating system for which they are usually installed by default. Given how fundamental C and Fortran are to scientific computing, it is unfortunate that Mac and Windows do not provide these compilers by default.
Detailed instructions for installing pomp and other software that we will use with it are provided in the following places:
Additional instructions on our course website
Please submit your solutions to Canvas as an Rmarkdown (.Rmd) file which the GSIs will compile into an HTML document. Your Rmd file can read in the Parus major data from the internet, e.g., by
dat <- read.csv("https://ionides.github.io/531w20/10/parus.csv")
Note: you will be using pomp version 2.x. There are some differences from pomp 1.x that may require attention if and when you look on the internet for pomp code.
Question 6.1. Reformulating the Ricker model.
The Ricker equation can be reparameterized so that the scaling of \(P_n\) is explicit: \[
P_{n+1} = r\,P_{n}\,\exp\left(-\frac{P_{n}}{k}\right).
\] Modify the pomp
object created in the notes to reflect this reparameterization. Also, Modify the measurement model so that the data \({y_n^*}\) is modeled as \[
Y_n |P_n \sim \mathrm{Negbin}(\phi\,P_n,\psi).
\] Here, \(\mathrm{Negbin}(\mu,\psi)\) is the negative binomial distribution with mean \(\mu\) and probability parameter \(\psi\), and therefore variance \(\mu/\psi\). This parameterization corresponds in R to rbinom(...,mu,prob)
. See ?rnbinom
for documentation on the negative binomial distribution and the R Extensions Manual section on distribution functions for information on how to access these in C.
Try simulating from a few choices of the parameters, and present one simulation from a set of parameters that shows oscillatory behavior.
Question 6.2. Coding a new POMP model.
Construct a pomp object for the Parus major data modeled using the stochastic Beverton-Holt model, \[ P_{n+1} = \frac{a\,P_n}{1+b\,P_n}\,\varepsilon_n, \] where \(a\) and \(b\) are parameters and \[ \varepsilon_t \sim \mathrm{Lognormal}(-\tfrac{1}{2}\sigma^2,\sigma^2). \] Assume the same measurement model as we used for the Ricker model. Try simulating from a few choices of the parameters. What are the similarities and differences between simulations you obtain from the Beverton-Holt model and those from the Ricker model? Present one simulation to support your comments.
Question 6.3. This feedback question is worth credit.
Explain which parts of your responses above made use of a source, meaning anything or anyone you consulted (including classmates or office hours) to help you write or check your answers. All sources are permitted. To encourage responsible use of these sources while maintaining class integrity, we require a response to this question, even if this may occasionally just say that you worked out everything entirely by yourself. See the syllabus for additional information on grading.
How long did this homework take? Report on any technical difficulties that arose.
This homework is conceptually quite simple, but involves overcoming various technical hurdles. The hurdles may be overcome quite quickly, or could turn into a longer battle.
To make progress on statistical inference for POMP models, we have to solve these underlying computational issues.
If you get stuck, ask for help from your peers and/or the GSIs and/or me. Please report how much time this homework ends up taking, to help me monitor how many difficulties are encountered.
The questions derive from material in a short course on Simulation-based Inference for Epidemiological Dynamics
Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.