RcppSMC GSoC Project

This page describes the progress made on the 2017 Google Summer of Code project RcppSMC, which is based on developing and enhancing the existing CRAN package RcppSMC.

Student:

Leah South

Mentors:

Adam Johansen and Dirk Eddelbuettel

Host organisation:

The R Project for Statistical Computing

Acknowledgements:

Thank you Adam and Dirk for giving me the chance to work on this exciting project. I've learnt a lot from both of you, especially on open source projects, C++ and particle filtering methods. I really appreciate the time you've spent explaining new concepts to me and discussing implementation details with me. I would also like to thank my PhD supervisor Chris Drovandi for encouraging my involvement in GSoC and for providing additional insights.

Background

Sequential Monte Carlo (SMC) methods use importance sampling and resampling to sample from a sequence of distributions. A common application of SMC is in signal processing, where SMC algorithms known as particle filters can be used to solve filtering problems. More recently, SMC methods have been used to perform inference in the context of static models. SMC can represent a powerful alternative to Markov chain Monte Carlo (MCMC) in this context because it is naturally adaptive, capable of sampling from complex and multimodal targets, and provides an estimate of the normalizing constant as a byproduct of the method.

Johansen (2009) developed the SMC template class (SMCTC) library in C++ for easy, efficient implementation of SMC methods. This library offers a single framework for SMC algorithms including particle filters and the more complex SMC samplers of Del Moral, Doucet and Jasra (2006). In 2012, Dirk Eddelbuettel and Adam Johansen developed RcppSMC, which provides access to SMCTC through R.

Prior to GSoC, RcppSMC provided support for a wide range of SMC methods and included three particle filtering examples. However, there were some commonly used and desirable SMC methods which were either not possible to implement with the existing library or required significant work by the user. The work done during GSoC was based around improving functionality and making common tasks easier for the user.

Progress

Details on the progress made during GSoC, including links to issue tickets and pull requests, are given below. These changes are described in more detail in the changelog. The differences between the code base prior to GSoC and at the completion of GSoC can be found here.

Work on the library:

• The packaging and interface were updated to use Rcpp attributes and the automatic RNGscope (4 and 5). This led to a simpler interface from R and simplified use of random number generators in the examples.
• Provided a generally applicable estimator for the ratio of normalizing constants and improved the existing estimators which are relevant to a small subset of cases (7). The normalizing constant is useful for Bayesian model choice and for pseudo-marginal methods where the likelihood can be estimated as the normalizing constant from a particle filter.
• Provided a framework for additional algorithm parameters and their adaptation (16, 19, 23 and 24). Being able to adapt the parameters based on the population of particles is a major benefit of SMC over standard MCMC, and this was not possible to do using the library prior to GSoC. Users can now adapt the parameters at several stages (including before an MCMC step).
• Added common adaptation methods for static Bayesian models (17, 20 and 21), including adaptation of the number of MCMC repeats, the covariance matrix for MCMC proposals and the temperature schedule.
• Switched to a population level object rather than a particle level object (2 and 3). This change was done early on as it affected many parts of the library. Some of the internal code was vectorized using RcppArmadillo, and further work on vectorization is planned (see future work section) to take advantage of the parallel nature of SMC.

Work on the examples:

• Added an example of SMC for static Bayesian models. Simple data annealing and likelihood annealing approaches to this example were both implemented (9 and 10). To demonstrate the new features for adaptation, which is essential for efficient implementation with static models, an adaptive version of likelihood annealing SMC was also implemented (18 and 22).
• Added an example to demonstrate how particle marginal Metropolis Hastings (and exact-approximate methods which use particle filtering estimates of the likelihood more generally) can be performed with the library (11 and 13). This example could not have been implemented with the existing library.

The repository connected to this wiki shows the progress on this work before formal pull requests were made.

Future Work

Some further developments which are planned for this package are given below.

• Multidimensional linear Gaussian models are often used when testing new methods because they can be solved exactly. Providing functions that return the exact Kalman filtering results and efficient approximations from the fully adapted auxiliary particle filter may be helpful.
• Generalize the examples. For example, change the linear regression example so that users can specify their own prior hyper-parameters, response vector and design matrix.
• Further improve the documentation by adding vignettes and potentially writing a journal article for JSS.
• There are many applications for which the particle values can be represented by a vector of real numbers. Future work may consider vectorizing calculations and allowing for whole population move and initialization functions to speed up this common use case.
• Add the option to track ancestral lines, which is useful for particle Gibbs (Andrieu, Doucet & Holenstein, 2010; Lindsten, Jordan and Schön, 2014).
• Add an implementation of the iterated auxiliary particle filter (iAPF, Guarniero, Johansen & Lee (2016)). The iAPF is more general than the auxiliary particle filter and implementing it will showcase how state of the art particle filtering methods can be implemented using RcppSMC.
• Potentially switch to Roxygen documentation for the .Rd files (8).

These plans will also be reflected in the package TODO file.

References

• Andrieu, C., Doucet, A., & Holenstein, R. (2010). Particle Markov chain Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3), 269-342.
• Del Moral, P., Doucet, A., & Jasra, A. (2006). Sequential Monte Carlo samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3), 411-436.
• Guarniero, P., Johansen, A. M., & Lee, A. (2016). The iterated auxiliary particle filter. Journal of the American Statistical Association (accepted).
• Johansen, A. M. (2009). SMCTC: sequential Monte Carlo in C++. Journal of Statistical Software, 30(6), 1-41.
• Lindsten, F., Jordan, M. I., & Schön, T. B. (2014). Particle Gibbs with ancestor sampling. Journal of Machine Learning Research, 15(1), 2145-2184.