PaperlayerView exact revision on arXiv

Source-linked AI summary

Statistical Inference for Partially Observed Markov Processes via the R Package pomp

Aaron A. King, Dao Nguyen, Edward L. Ionides

arXiv:1509.00503v2stat.ME

TL;DR

POMP models are widely used for noisy observations of latent Markov processes, yet broad nonlinear and non-Gaussian models require specialized computational methods. The paper presents pomp as a general R framework implementing diverse Monte Carlo methods and demonstrates its use across toy and epidemiological models. The package supports model-independent algorithm development and both Bayesian and non-Bayesian analyses, while examples show that method performance and comparisons can depend on settings and problem context.

  • Problem

    Broad classes of nonlinear, non-Gaussian POMP models need flexible software that can support the growing range of specialized statistical methodologies.

  • Method

    The paper presents pomp as a general, object-oriented R framework separating POMP model specification from inference algorithms and implementing methods including SMC, PMCMC, ABC, synthetic likelihood, forecasting, and trajectory matching.

  • Results

    pomp makes a wide range of inference algorithms available from one model implementation and supports objective comparison of alternative models and methods.

  • Takeaways & Limitations

    pomp serves both as a POMP data-analysis tool and as a platform for developing and extending inference algorithms, including plug-and-play methods.

  • Takeaways & Limitations

    The examples are limited in generalizability, and method comparisons can depend on algorithmic settings and the problem of interest.

Abstract

from arXiv · show

Partially observed Markov process (POMP) models, also known as hidden Markov models or state space models, are ubiquitous tools for time series analysis. The R package pomp provides a very flexible framework for Monte Carlo statistical investigations using nonlinear, non-Gaussian POMP models. A range of modern statistical methods for POMP models have been implemented in this framework including sequential Monte Carlo, iterated filtering, particle Markov chain Monte Carlo, approximate Bayesian computation, maximum synthetic likelihood estimation, nonlinear forecasting, and trajectory matching. In this paper, we demonstrate the application of these methodologies using some simple toy problems. We also illustrate the specification of more complex POMP models, using a nonlinear epidemiological model with a discrete population, seasonality, and extra-demographic stochasticity. We discuss the specification of user-defined models and the development of additional methods within the programming environment provided by pomp.

1. Introduction

POMP models describe noisy, incomplete observations of latent Markov processes, but software has struggled to support their broad model classes and statistical methods. pomp addresses this by providing a general framework focused on nonlinear, non-Gaussian models with large state spaces.

  • POMP models represent incomplete and noisy measurements of a latent, unobserved Markov process.
  • pomp uses a general and abstract representation so algorithms implemented within it apply to arbitrary POMP models.
  • The package focuses on nonlinear, non-Gaussian POMP models with large state spaces.
  • Methods exploiting linear-Gaussian approximations, small discrete state spaces, or exact hidden-Markov-model structure are not yet part of pomp.
  • pomp supports both Bayesian and non-Bayesian data analyses by allowing parameter priors and implementing several Bayesian methods.

2. POMP models and their representation in pomp

A POMP model is defined through latent-process and measurement components and represented in pomp as an object whose slots encode those components. This representation supports simulation, inference, initialization choices, and models that vary with time or observed covariates.

  • A POMP model combines an initial density, conditional transition density, and measurement density to determine the joint process density.
  • The transition and measurement densities may depend explicitly on the observation index n, allowing nonhomogeneous POMP models.
  • pomp represents each model as an S4 object of class ‘pomp’, with slots encoding model components used by class methods.
  • The rprocess, dprocess, rmeasure, and dmeasure arguments specify transition probabilities and measurement densities, with dprocess unnecessary for plug-and-play methods.
  • Initial latent states can be fixed by experimental design, set to a stationary distribution, or estimated as parameters when no scientifically motivated initialization exists.
  • Observed covariate processes can enter the arbitrary model densities through the covar argument, supporting conditional inference for nonhomogeneous POMPs.

3. Methodology for POMP models

The pomp framework supports inference for broad POMP models by combining model representations with Monte Carlo methods, especially plug-and-play algorithms. Its methodology spans likelihood-based, Bayesian, and feature-based approaches, including particle filtering, iterated filtering, and PMCMC.

  • Methodological framework: POMP models are analyzed through full-information or feature-based methods, and through Bayesian or frequentist, plug-and-play or non-plug-and-play categories.Feature-based methods emphasize selected data features rather than the full likelihood and may lose statistical efficiency when low-dimensional sufficient statistics are unavailable.
  • Methodological framework: Pomp focuses on nonlinear, non-Gaussian POMP models with large state spaces, where methods leveraging POMP structure are needed.The package is flexible enough to encode arbitrary POMP models and methods, but its development has focused on plug-and-play methods.
  • Sequential Monte Carlo: Sequential Monte Carlo propagates particles with the dynamic model, assimilates observations, and estimates the likelihood through a particle-filter representation of the latent process.The basic particle filter is implemented through pfilter and has O(J) complexity per observation step in the supplied algorithm description.
  • Sequential Monte Carlo: Basic SMC can suffer particle depletion when observations are extremely unlikely under the model, while systematic resampling reduces Monte Carlo variability relative to multinomial resampling.The paper notes that investigating alternative models may be preferable to relying only on algorithmic elaborations to address severe depletion.
  • Iterated filtering: Iterated filtering maximizes an SMC-based likelihood by temporarily allowing parameters to follow a random walk whose intensity decreases across filtering repetitions.The added variability smooths the likelihood surface, combats particle depletion, supports gradient estimation, and preserves plug-and-play operation.
  • Particle Markov chain Monte Carlo: Particle Markov chain Monte Carlo combines SMC likelihood evaluation with parameter-space MCMC moves to enable full-information plug-and-play Bayesian inference.Particle marginal Metropolis-Hastings uses the unbiased SMC likelihood estimate within the Metropolis-Hastings update to target the desired posterior.

4. Model construction and data analysis: Simple examples

The section builds a Gompertz POMP model in pomp by specifying process and measurement components, then simulates and visualizes data. The model’s log transformation also permits exact likelihood calculations with the Kalman filter for comparison with general methods.

  • Gompertz model: The Gompertz example specifies a discrete-time population process with lognormal process noise and lognormal measurement error.K is carrying capacity, r is positive, and process noise has log ε_t ∼ Normal(0, σ^2).
  • Likelihood comparison: After logarithmic transformation, the Gompertz model is linear and Gaussian, so the Kalman filter supplies exact likelihood values for comparison.The example therefore contrasts generally applicable plug-and-play methods with an exact method available for this special model.
  • Model construction: A pomp object encodes data, observation times, a process simulator, and an observation simulator through the pomp constructor.The process simulator advances the latent state by one time step, while the measurement simulator generates an observation from the observation process.
  • Model construction: The measurement density evaluates the likelihood or log likelihood of an observation given the latent state and parameters, enabling pfilter, mif, and pmcmc.The log argument selects whether the density or log density is returned.
  • Simulation and visualization: The model can be simulated and plotted, replacing missing observations with simulated data and saving the simulation parameters internally.Figure 1 shows the resulting simulated observations from the Gompertz model.

4.2. Computing likelihood using SMC

This section uses sequential Monte Carlo through pfilter to estimate POMP likelihoods from particle simulations. In the Gompertz example, the particle estimate is close to the exact Kalman-filter likelihood, with discrepancy attributed to Monte Carlo error.

  • Sequential Monte Carlo: The particle filter pfilter estimates a POMP likelihood using a chosen number of particles, trading lower Monte Carlo error against greater computational burden.The implementation requires process simulation and measurement-density evaluation for this full-information plug-and-play method.
  • Gompertz example: With 1000 particles, pfilter estimates the Gompertz log likelihood at the true parameters as 36.27102.The estimate is obtained by applying pfilter and then logLik to the resulting object.
  • Gompertz example: At a parameter guess with r, K, and σ each 50% above their true values, the estimated log likelihood is 25.19585.The comparison evaluates the likelihood at a deliberately displaced point in parameter space.
  • Gompertz example: 36.01 is the exact Kalman-filter log likelihood at the true parameters, compared with 36.27 from the 1000-particle filter.The difference is attributed to Monte Carlo error because the particle filter provides an unbiased likelihood estimate.
  • Monte Carlo error: Repeating pfilter and averaging likelihood estimates can reduce Monte Carlo error while also permitting estimation of that error.The averaging is performed on likelihoods rather than log likelihoods because the particle estimate is unbiased on the likelihood scale.

4.3. Maximum likelihood estimation via iterated filtering

This section applies iterated filtering to estimate Gompertz parameters using repeated particle-filter searches from multiple starting points. The resulting estimate is close to the exact MLE, with likelihood agreement within about 0.1 log units.

  • Parameter estimation: Iterated filtering estimates positive Gompertz parameters on an unconstrained transformed scale using logarithmic transformations and their inverses.The transformations enforce positivity during estimation.
  • Parameter estimation: The analysis runs 10 iterated-filtering searches from dispersed starting points and estimates r, σ, and τ.The computations are parallelized with foreach, and likelihoods are repeatedly evaluated with particle filtering.
  • Convergence and selection: The 10 mif runs converge to different point estimates, so the analysis focuses on the run with the highest estimated likelihood.Likelihoods are evaluated repeatedly to reduce Monte Carlo error before selecting the estimate.
  • Convergence and selection: Particle-filter likelihoods are averaged rather than log likelihoods because the particle filter produces an unbiased likelihood estimate.Repeated evaluations are used to reduce Monte Carlo error in the likelihood comparison.
  • Results: Table 3 compares mif estimates with the exact MLE and true parameter values, reporting SMC likelihood estimates, standard errors, and the exact likelihood.Both the mif MLE and the truth lie within the ideal likelihood-ratio 95% confidence set, and the mif MLE is close to the exact MLE.
  • Results: The mif procedure maximizes the likelihood up to an error of about 0.1 log units in the Gompertz example.This comparison uses the exact Kalman-filter likelihood available for the tractable transformed model.

4.4. Full-information Bayesian inference via PMCMC

The example specifies a prior density for Bayesian inference and uses PMCMC with five independent random-walk chains after locating an MLE neighborhood. Diagnostics show limited effective sample size for r, while PMCMC is computationally less efficient than iterated filtering.

  • Bayesian setup: The Gompertz example specifies dprior for random-walk Metropolis-Hastings, which requires prior-density evaluation but not prior simulation.The prior density is implemented using uniform distributions bounded by hyperparameters.
  • PMCMC analysis: PMCMC draws posterior samples using a multivariate normal random-walk proposal after iterated filtering locates a neighborhood of the MLE.The analysis runs five independent chains with diagonal proposal variance.
  • Diagnostics: 250 is the lowest combined effective sample size, for r, after 40000 proposal steps.The diagnostic figure reports this as a modest number of proposal steps for that parameter.
  • Computational considerations: PMCMC can require orders of magnitude more computation than iterated filtering and cannot be reliably assessed from a single chain.The paper recommends multiple runs from a range of starting points and uses diagnostic plots to check convergence.

4.5. A second example: The Ricker model

The Ricker example introduces a discrete-time stochastic population model with noisy observations and demonstrates how to implement it in pomp. The model uses compiled C snippets and parameter transformations within a class ‘pomp’ object.

  • Ricker model: The Ricker model describes population size N_t in a discrete-time population model and assumes that measurements Y_t are noisy.The measurement process is described separately from the latent population process.
  • Model implementation: pomp supports C snippets for the Ricker state-process simulator, measurement simulator, and measurement density, which can provide manyfold speed-ups over R.The snippets use stochastic population evolution and Poisson measurement noise.
  • Model implementation: The implementation treats both N and e as state variables, while pomp automatically constructs variable declarations.The give_log flag selects likelihood or log-likelihood output for the measurement density.
  • Parameter handling: Logarithmic and exponential transformations are added to enforce parameter constraints in the Ricker model.The transformations cover r, sigma, phi, and N.0.
  • Model construction: The complete Ricker model object combines discrete-time simulation, measurement functions, parameter transformations, parameter names, state names, and observation times.The object is constructed with pomp and initialized with values for r, sigma, phi, N.0, and e.0.

4.6. Feature-based synthetic likelihood maximization

The section demonstrates feature-based synthetic likelihood maximization for the Ricker model using probes applied to observed and simulated data. Comparison with iterated filtering shows that full-likelihood maximization is tractable, whereas synthetic likelihood is statistically less efficient.

  • Probes: pomp probes are functions that map real or simulated data to scalar or vector quantities, supporting feature-based comparisons.The example uses marginal, autocorrelation, and nonlinear autoregressive probes.
  • Probe diagnostics: The probe workflow applies the probe functions to observed data and 1000 simulated data sets, then produces summaries and diagnostic plots.The plots assess model-data agreement and correlations among probes.
  • Synthetic likelihood optimization: Probe matching uses optimization to identify parameters that maximize the synthetic likelihood returned by the probe analysis.The example uses Nelder-Mead optimization for r, sigma, and phi.
  • Results: Likelihood maximization for the Ricker model is possible and not difficult, contrary to Wood’s claim that the full likelihood was intractable.Table 4 compares the guess, truth, MLE, and MSLE using Monte Carlo likelihoods and synthetic likelihoods.
  • Results: Maximum synthetic likelihood is statistically inefficient relative to likelihood because synthetic likelihood discards some information in the data.At the MSLE, the estimated log likelihood is smaller than at the truth.

4.7. Bayesian feature matching via ABC

The ABC example uses a single simulation per iteration, making iterations quicker than synthetic likelihood evaluations but requiring feature scaling in advance. In this example ABC mixes somewhat faster than PMCMC, while producing broader posterior distributions and potentially losing information.

  • ABC computation: ABC uses one simulation per iteration, so each iteration is essentially the cost of SMC with one particle or synthetic likelihood with one simulation.Because ABC cannot determine relative feature scaling during each evaluation, scaling must be supplied in advance.
  • ABC setup: The example runs ABC with probes, a tolerance of 2, precomputed scaling, and 4e6 iterations across five chains.The probe set includes mean, autocorrelation, and marginal-distribution features.
  • Results: 450 is the lowest effective sample size for ABC’s r parameter, compared with 250 for PMCMC under matched total computational effort.The paper concludes that ABC mixes somewhat more rapidly in this example.
  • Results: ABC produces somewhat broader posterior distributions than the full-information PMCMC posteriors in this example.The comparison concerns the statistical efficiency of ABC on this model and data set.
  • Scope and caveats: The conclusions may depend on algorithmic settings and a particular model-data example, while ABC may remain the only practical Bayesian method in some situations.The paper warns that substantial information loss can occur even with scientifically reasoned features.

4.8. Parameter estimation by simulated quasi-likelihood

The paper compares nonlinear forecasting with iterated filtering for simulated quasi-likelihood estimation, finding that nonlinear forecasting is faster but can perform poorly under likelihood-based evaluation.

  • A simulation study compares nlf with mif on correctly specified Gompertz-model data, while noting that it does not address fitting misspecified models.
  • Figure 6 compares mif’s maximum likelihood estimate with nlf’s maximum simulated quasi-likelihood estimate across 10 simulated datasets.Panel A evaluates estimated log likelihood, while panel B evaluates simulated log quasi-likelihood.
  • The MSQL estimate can fall many units of log likelihood short of the MLE, although likelihood-based inference is nearly as good at optimizing nlf’s target criterion.
  • Each mif optimization took 26.9 sec versus 3.7 sec for nlf, but extra computation or algorithmic adjustments might change either method’s performance.

5. A more complex example: Epidemics in continuous time

The paper develops a continuous-time stochastic SIR POMP model with discrete demographic events and noisy case reporting, then implements it in pomp using simulation and measurement components.

  • The SIR model divides hosts into susceptible, infected, and recovered classes, with births, deaths, infection, and recovery changing class sizes through time.
  • The force of infection is λ(t) = β I(t)/P, while equal birth and death rates imply a constant expected population size.
  • The deterministic skeleton is obtained by taking expectations of the continuous-time process equations and letting the time increment approach zero.
  • Reported cases are modeled from accumulated true incidence using a negative binomial measurement process with reporting rate ρ and size parameter θ.
  • The implementation uses pomp Csnippets for model components, enforces integer initial states summing to population size, and produces simulated data and latent states.The constructed models use rprocess, rmeasure, and dmeasure components; Fig. 9 displays one resulting simulation.
  • The approximate tau-leap implementation holds transition rates constant over an Euler interval and uses multinomial transitions for individuals leaving each class.

6. Conclusion

The conclusion presents pomp as both a general POMP data-analysis tool and an extensible platform for inference-algorithm development, with support for computationally demanding applications.

  • pomp separates model specification from inference algorithms, making alternative models and methods easier to compare within one framework.
  • Its general model specification language makes a wide range of inference algorithms available for POMP analyses.
  • The platform is especially convenient for plug-and-play methods while also accommodating algorithms using dprocess and deterministic skeleton components.
  • Complex models and large datasets can require C-based acceleration and multiprocessor computing for computationally intensive tasks.
  • The package supplies examples, interactive demonstrations, documentation, and an introductory tutorial to support new model implementations.
Loading 1509.00503v2…