Pierre E. Jacob

Using parallel computation to improve Independent Metropolis-Hastings based estimation

In this article, we consider the implications of the fact that parallel raw-power can be exploited by a generic Metropolis–Hastings algorithm if the proposed values are independent from the current value of the Markov chain. In particular, we present improvements to the independent Metropolis–Hastings algorithm that significantly decrease the variance of any estimator derived from the MCMC output, at a null computing cost since those improvements are based on a fixed number of target density evaluations that can be produced in parallel. The techniques developed in this article do not jeopardize the Markovian convergence properties of the algorithm, since they are based on the Rao–Blackwell principles of Gelfand and Smith (1990), already exploited in the work of Casella and Robert (1996), Atchadé and Perron (2005), and Douc and Robert (2011). We illustrate those improvements both on a toy normal example and on a classical probit regression model, but stress the fact that they are applicable in any case where the independent Metropolis–Hastings is applicable. The code used in this article is available as supplementary material.

Parallel Resampling in the Particle Filter

Modern parallel computing devices, such as the graphics processing unit (GPU), have gained significant traction in scientific and statistical computing. They are particularly well-suited to data-parallel algorithms such as the particle filter, or more generally sequential Monte Carlo (SMC), which are increasingly used in statistical inference. SMC methods carry a set of weighted particles through repeated propagation, weighting, and resampling steps. The propagation and weighting steps are straightforward to parallelize, as they require only independent operations on each particle. The resampling step is more difficult, as standard schemes require a collective operation, such as a sum, across particle weights. Focusing on this resampling step, we analyze two alternative schemes that do not involve a collective operation (Metropolis and rejection resamplers), and compare them to standard schemes (multinomial, stratified, and systematic resamplers). We find that, in certain circumstances, the alternative resamplers can perform significantly faster on a GPU, and to a lesser extent on a CPU, than the standard approaches. Moreover, in single precision, the standard approaches are numerically biased for upward of hundreds of thousands of particles, while the alternatives are not. This is particularly important given greater single- than double-precision throughput on modern devices, and the consequent temptation to use single precision with a greater number of particles. Finally, we provide auxiliary functions useful for implementation, such as for the permutation of ancestry vectors to enable in-place propagation. Supplementary materials are available online.

Path storage in the particle filter

This article considers the problem of storing the paths generated by a particle filter and more generally by a sequential Monte Carlo algorithm. It provides a theoretical result bounding the expected memory cost by T+CNlogN where T is the time horizon, N is the number of particles and C is a constant, as well as an efficient algorithm to realise this. The theoretical result and the algorithm are illustrated with numerical experiments.

Feynman-Kac particle integration with geometric interacting jumps

This article is concerned with the design and analysis of discrete time Feynman-Kac particle integration models with geometric interacting jump processes. We analyze two general types of model, corresponding to whether the reference process is in continuous or discrete time. For the former, we consider discrete generation particle models defined by arbitrarily fine time mesh approximations of the Feynman-Kac models with continuous time path integrals. For the latter, we assume that the discrete process is observed at integer times and we design new approximation models with geometric interacting jumps in terms of a sequence of intermediate time steps between the integers. In both situations, we provide nonasymptotic bias and variance theorems w.r.t. the time step and the size of the system, yielding what appear to be the first results of this type for this class of Feynman-Kac particle integration models. We also discuss uniform convergence estimates w.r.t. the time horizon. Our approach is based on an original semigroup analysis with first order decompositions of the fluctuation errors.

Smoothing with Couplings of Conditional Particle Filters

In state space models, smoothing refers to the task of estimating a latent stochastic process given noisy measurements related to the process. We propose an unbiased estimator of smoothing expectat ...

Bayesian Inference in Non-Markovian State-Space Models With Applications to Battery Fractional-Order Systems

Battery impedance spectroscopy models are given by fractional-order (FO) differential equations. In the discrete-time domain, they give rise to state-space models where the latent process is not Markovian. Parameter estimation for these models is, therefore, challenging, especially for noncommensurate FO models. In this paper, we propose a Bayesian approach to identify the parameters of generic FO systems. The computational challenge is tackled with particle Markov chain Monte Carlo methods, with an implementation specifically designed for the non-Markovian setting. Two examples are provided. In a first example, the approach is applied to identify a battery commensurate FO model with a single constant phase element (CPE) by using real data. We compare the proposed approach to an instrumental variable method. Then, we consider a noncommensurate FO model with more than one CPE and synthetic data sets, investigating how the proposed method enables the study of various effects on parameter identification, such as the data length, the magnitude of the input signal, the choice of prior, and the measurement noise.

Bayesian model comparison with the Hyvärinen score: computation and consistency

ABSTRACT The Bayes factor is a widely used criterion in model comparison and its logarithm is a difference of out-of-sample predictive scores under the logarithmic scoring rule. However, when some of the candidate models involve vague priors on their parameters, the log-Bayes factor features an arbitrary additive constant that hinders its interpretation. As an alternative, we consider model comparison using the Hyvärinen score. We propose a method to consistently estimate this score for parametric models, using sequential Monte Carlo methods. We show that this score can be estimated for models with tractable likelihoods as well as nonlinear non-Gaussian state-space models with intractable likelihoods. We prove the asymptotic consistency of this new model selection criterion under strong regularity assumptions in the case of nonnested models, and we provide qualitative insights for the nested case. We also use existing characterizations of proper scoring rules on discrete spaces to extend the Hyvärinen score to discrete observations. Our numerical illustrations include Lévy-driven stochastic volatility models and diffusion models for population dynamics. Supplementary materials for this article are available online.