Sequential Monte Carlo Sampling Research Articles

AstroABC : An Approximate Bayesian Computation Sequential Monte Carlo sampler for cosmological parameter estimation

Given the complexity of modern cosmological parameter inference where we are faced with non-Gaussian data and noise, correlated systematics and multi-probe correlated data sets, the Approximate Bayesian Computation (ABC) method is a promising alternative to traditional Markov Chain Monte Carlo approaches in the case where the Likelihood is intractable or unknown. The ABC method is called "Likelihood free" as it avoids explicit evaluation of the Likelihood by using a forward model simulation of the data which can include systematics. We introduce astroABC, an open source ABC Sequential Monte Carlo (SMC) sampler for parameter estimation. A key challenge in astrophysics is the efficient use of large multi-probe datasets to constrain high dimensional, possibly correlated parameter spaces. With this in mind astroABC allows for massive parallelization using MPI, a framework that handles spawning of jobs across multiple nodes. A key new feature of astroABC is the ability to create MPI groups with different communicators, one for the sampler and several others for the forward model simulation, which speeds up sampling time considerably. For smaller jobs the Python multiprocessing option is also available. Other key features include: a Sequential Monte Carlo sampler, a method for iteratively adapting tolerance levels, local covariance estimate using scikit-learn's KDTree, modules for specifying optimal covariance matrix for a component-wise or multivariate normal perturbation kernel, output and restart files are backed up every iteration, user defined metric and simulation methods, a module for specifying heterogeneous parameter priors including non-standard prior PDFs, a module for specifying a constant, linear, log or exponential tolerance level, well-documented examples and sample scripts. This code is hosted online at https://github.com/EliseJ/astroABC

Inferring Volatility Dynamics Using Stock Prices and Variance Swap Rates: An Efficient Bayesian Approach

We propose a particle-based Bayesian method based on a quasi-Monte-Carlo particle filter and a two-stage sequential Monte Carlo sampler with common random numbers to estimate a class of volatility models that take into account co-jumps of underlying prices and volatility, non-affineness, and self-excitation. This method is robust to the presence of outliers in the data and can deliver almost full-likelihood-based inference at a much smaller computing cost. Empirical study using both S&P 500 index and variance swap rates shows that both non-affineness and self-excitation need to be introduced in modeling volatility dynamics.

A Bayesian approach to multiscale inverse problems with on-the-fly scale determination

A Bayesian computational approach is presented to provide a multi-resolution estimate of an unknown spatially varying parameter from indirect measurement data. In particular, we are interested in spatially varying parameters with multiscale characteristics. In our work, we consider the challenge of not knowing the characteristic length scale(s) of the unknown a priori, and present an algorithm for on-the-fly scale determination. Our approach is based on representing the spatial field with a wavelet expansion. Wavelet basis functions are hierarchically structured, localized in both spatial and frequency domains and tend to provide sparse representations in that a large number of wavelet coefficients are approximately zero. For these reasons, wavelet bases are suitable for representing permeability fields with non-trivial correlation structures. Moreover, the intra-scale correlations between wavelet coefficients form a quadtree, and this structure is exploited to identify additional basis functions to refine the model. Bayesian inference is performed using a sequential Monte Carlo (SMC) sampler with a Markov Chain Monte Carlo (MCMC) transition kernel. The SMC sampler is used to move between posterior densities defined on different scales, thereby providing a computationally efficient method for adaptive refinement of the wavelet representation. We gain insight from the marginal likelihoods, by computing Bayes factors, for model comparison and model selection. The marginal likelihoods provide a termination criterion for our scale determination algorithm. The Bayesian computational approach is rather general and applicable to several inverse problems concerning the estimation of a spatially varying parameter. The approach is demonstrated with permeability estimation for groundwater flow using pressure sensor measurements.

Multilevel sequential Monte Carlo samplers

In this article we consider the approximation of expectations w.r.t. probability distributions associated to the solution of partial differential equations (PDEs); this scenario appears routinely in Bayesian inverse problems. In practice, one often has to solve the associated PDE numerically, using, for instance finite element methods which depend on the step-size level hL. In addition, the expectation cannot be computed analytically and one often resorts to Monte Carlo methods. In the context of this problem, it is known that the introduction of the multilevel Monte Carlo (MLMC) method can reduce the amount of computational effort to estimate expectations, for a given level of error. This is achieved via a telescoping identity associated to a Monte Carlo approximation of a sequence of probability distributions with discretization levels ∞>h0>h1⋯>hL. In many practical problems of interest, one cannot achieve an i.i.d. sampling of the associated sequence and a sequential Monte Carlo (SMC) version of the MLMC method is introduced to deal with this problem. It is shown that under appropriate assumptions, the attractive property of a reduction of the amount of computational effort to estimate expectations, for a given level of error, can be maintained within the SMC context. That is, relative to exact sampling and Monte Carlo for the distribution at the finest level hL. The approach is numerically illustrated on a Bayesian inverse problem.

Fast Sampling for Time-Varying Determinantal Point Processes

Determinantal Point Processes (DPPs) are stochastic models which assign each subset of a base dataset with a probability proportional to the subset’s degree of diversity. It has been shown that DPPs are particularly appropriate in data subset selection and summarization (e.g., news display, video summarizations). DPPs prefer diverse subsets while other conventional models cannot offer. However, DPPs inference algorithms have a polynomial time complexity which makes it difficult to handle large and time-varying datasets, especially when real-time processing is required. To address this limitation, we developed a fast sampling algorithm for DPPs which takes advantage of the nature of some time-varying data (e.g., news corpora updating, communication network evolving), where the data changes between time stamps are relatively small. The proposed algorithm is built upon the simplification of marginal density functions over successive time stamps and the sequential Monte Carlo (SMC) sampling technique. Evaluations on both a real-world news dataset and the Enron Corpus confirm the efficiency of the proposed algorithm.

Toward Automatic Model Comparison: An Adaptive Sequential Monte Carlo Approach

Model comparison for the purposes of selection, averaging, and validation is a problem found throughout statistics. Within the Bayesian paradigm, these problems all require the calculation of the posterior probabilities of models within a particular class. Substantial progress has been made in recent years, but difficulties remain in the implementation of existing schemes. This article presents adaptive sequential Monte Carlo (SMC) sampling strategies to characterize the posterior distribution of a collection of models, as well as the parameters of those models. Both a simple product estimator and a combination of SMC and a path sampling estimator are considered and existing theoretical results are extended to include the path sampling variant. A novel approach to the automatic specification of distributions within SMC algorithms is presented and shown to outperform the state of the art in this area. The performance of the proposed strategies is demonstrated via an extensive empirical study. Comparisons with state-of-the-art algorithms show that the proposed algorithms are always competitive, and often substantially superior to alternative techniques, at equal computational cost and considerably less application-specific implementation effort. Supplementary materials for this article are available online.

Sequential Monte Carlo Sampling for State Space Models

Sequential Monte Carlo Sampling for State Space Models

Striated Metropolis–Hastings sampler for high-dimensional models

Having efficient and accurate samplers for simulating the posterior distribution is crucial for Bayesian analysis. We develop a generic posterior simulator called the “dynamic striated Metropolis–Hastings (DSMH)” sampler. Grounded in the Metropolis–Hastings algorithm, it pools the strengths from the equi-energy and sequential Monte Carlo samplers while avoiding the weaknesses of the standard Metropolis–Hastings algorithm and those of importance sampling. In particular, the DSMH sampler possesses the capacity to cope with extremely irregular distributions that contain winding ridges and multiple peaks; and it is robust to how the sampling procedure progresses across stages. The high-dimensional application studied in this paper provides a natural platform for testing any generic sampler.

FORWARD AND INVERSE UNCERTAINTY QUANTIFICATION USING MULTILEVEL MONTE CARLO ALGORITHMS FOR AN ELLIPTIC NONLOCAL EQUATION

This paper considers uncertainty quantification for an elliptic nonlocal equation. In particular, it is assumed that the parameters which define the kernel in the nonlocal operator are uncertain and a priori distributed according to a probability measure. It is shown that the induced probability measure on some quantities of interest arising from functionals of the solution to the equation with random inputs is well-defined; as is the posterior distribution on parameters given observations. As the elliptic nonlocal equation cannot be solved approximate posteriors are constructed. The multilevel Monte Carlo (MLMC) and multilevel sequential Monte Carlo (MLSMC) sampling algorithms are used for a priori and a posteriori estimation, respectively, of quantities of interest. These algorithms reduce the amount of work to estimate posterior expectations, for a given level of error, relative to Monte Carlo and i.i.d. sampling from the posterior at a given level of approximation of the solution of the elliptic nonlocal equation.

TRANSITIONAL ANNEALED ADAPTIVE SLICE SAMPLING FOR GAUSSIAN PROCESS HYPER-PARAMETER ESTIMATION

Surrogate models have become ubiquitous in science and engineering for their capability of emulating expensive computer codes, necessary to model and investigate complex phenomena. Bayesian emulators based on Gaussian processes adequately quantify the uncertainty that results from the cost of the original simulator, and thus the inability to evaluate it on the whole input space. However, it is common in the literature that only a partial Bayesian analysis is carried out, whereby the underlying hyper-parameters are estimated via gradient-free optimisation or genetic algorithms, to name a few methods. On the other hand, maximum a posteriori (MAP) estimation could discard important regions of the hyper-parameter space. In this paper, we carry out a more complete Bayesian inference, that combines Slice Sampling with some recently developed Sequential Monte Carlo samplers. The resulting algorithm improves the mixing in the sampling through delayed-rejection, the inclusion of an annealing scheme akin to Asymptotically Independent Markov Sampling and parallelisation via Transitional Markov Chain Monte Carlo. Examples related to the estimation of Gaussian process hyper-parameters are presented. For the purpose of reproducibility, further development, and use in other applications, the code to generate the examples in this paper is freely available for download at http://github.com/agarbuno/ta2s2_codes

A probabilistic framework for the exploration of enzymatic capabilities based on feasible kinetics and control analysis

BackgroundAnalysis of limiting steps within enzyme-catalyzed reactions is fundamental to understand their behavior and regulation. Methods capable of unravelling control properties and exploring kinetic capabilities of enzymatic reactions would be particularly useful for protein and metabolic engineering. While single-enzyme control analysis formalism has previously been applied to well-studied enzymatic mechanisms, broader application of this formalism is limited in practice by the limited amount of kinetic data and the difficulty of describing complex allosteric mechanisms. MethodsTo overcome these limitations, we present here a probabilistic framework enabling control analysis of previously unexplored mechanisms under uncertainty. By combining a thermodynamically consistent parameterization with an efficient Sequential Monte Carlo sampler embedded in a Bayesian setting, this framework yields insights into the capabilities of enzyme-catalyzed reactions with modest kinetic information, provided that the catalytic mechanism and a thermodynamic reference point are defined. ResultsThe framework was used to unravel the impact of thermodynamic affinity, substrate saturation levels and effector concentrations on the flux control and response coefficients of a diverse set of enzymatic reactions. ConclusionsOur results highlight the importance of the metabolic context in the control analysis of isolated enzymes as well as the use of statistically sound methods for their interpretation. General SignificanceThis framework significantly expands our current capabilities for unravelling the control properties of general reaction kinetics with limited amount of information. This framework will be useful for both theoreticians and experimentalists in the field.

Efficient Sequential Monte-Carlo Samplers for Bayesian Inference

In many problems, complex non-Gaussian and/or nonlinear models are required to accurately describe a physical system of interest. In such cases, Monte-Carlo algorithms are remarkably flexible and extremely powerful approaches to solve such inference problems. However, in the presence of a high-dimensional and/or multimodal posterior distribution, it is widely documented that standard Monte-Carlo techniques could lead to poor performance. In this paper, the study is focused on a Sequential Monte-Carlo (SMC) sampler framework, a more robust and efficient Monte-Carlo algorithm. Although this approach presents many advantages over traditional Monte-Carlo methods, the potential of this emergent technique is, however, largely underexploited in signal processing. In this paper, we aim at proposing some novel strategies to improve the efficiency and facilitate practical implementation of the SMC sampler. First, we propose an automatic and adaptive strategy that selects the sequence of distributions within the SMC sampler that minimizes the asymptotic variance of the estimator of the posterior normalization constant. The second original contribution we present improves the global efficiency of the SMC sampler by introducing a novel correction mechanism that allows the use of the particles generated through all of the iterations of the algorithm (instead of only particles from the last iteration). This is a significant contribution as it removes the need to discard a large portion of the samples obtained, as is standard in standard SMC methods. This will improve estimation performance in practical settings where the computational budget is important to consider.

A Bayesian Perspective on Multiple Source Localization in Wireless Sensor Networks

In this paper, we address the challenging problem of multiple source localization in wireless sensor networks (WSN). We develop an efficient statistical algorithm, based on the novel application of sequential Monte Carlo (SMC) sampler methodology, that is able to deal with an unknown number of sources given quantized data obtained at the fusion center from different sensors with imperfect wireless channels. We also derive the posterior Cramér-Rao bound (PCRB) of the source location estimate. The PCRB is used to analyze the accuracy of the proposed SMC sampler algorithm and the impact that quantization has on the accuracy of location estimates of the sources. Extensive experiments show the benefits of the proposed scheme in terms of the accuracy of the estimation method that is required for model selection (i.e., the number of sources) and the estimation of the source characteristics compared to the classical importance sampling method.

Iterative ensemble smoothers in the annealed importance sampling framework

Iterative ensemble techniques for solving inverse problems has recently gained a lot of interest in many geophysical communities. This popularity is attributed to the simplicity of implementation, general reliability and the ability to deal with the forward model as a black box without requiring the implementation of analytical gradients. Although several variants exist, we focus on the ensemble smoother with multiple data assimilation. This study highlights the similarity between the ensemble smoother and other existing techniques such as particle flow and annealed importance sampling. It is shown how a sequential Monte Carlo sampler can be used in combination with an annealing process to weight-correct the sampling procedure used in the ensemble smoother. Two different approximations in high dimensions, where the curse of dimensionality is unavoidable, are also presented. The methods proposed are compared with an MCMC run on a synthetic reservoir model.

Sequential Monte-Carlo sampling based on a committee of artificial neural networks for posterior state estimation and residual lifetime prediction

The application of Bayesian methods to the problem of fatigue crack growth prediction has been growing in recent years. In particular, sequential Monte-Carlo sampling is often presented as an efficient model-based technique to filter the sequential measures of the damage evolution provided as an input to the algorithm. However, a lot of measures are required to reliably identify the system state condition and the underlying model parameters. Many studies rely on the availability of a relatively dense sequence of crack length measures during damage evolution, made in most cases impractical by the consequent maintenance costs. Thus, real-time damage diagnosis is a requirement to enable prognostic health management.This work focuses on the application of sequential Monte-Carlo sampling to estimate the probabilistic residual life of a structural component subjected to fatigue crack propagation, while real-time estimation of crack length is provided through a committee of artificial neural networks, trained with finite element simulated strain patterns. Multiple crack length observations are available at each discrete time and are provided as the input to the prognostic system, based on a sequential importance resampling algorithm. Each time a new set of measures is available, the algorithm evaluates the posterior distribution of the augmented state vector, including the crack length and a material parameter governing damage evolution. This filtered information is used to numerically update the probability density functions of the residual life of the component. The methodology is applied first to a simulated crack and then to a metallic stiffened panel specimen subject to fatigue crack growth.

Real-Time Prognosis of Crack Growth Evolution Using Sequential Monte Carlo Methods and Statistical Model Parameters

A probabilistic method to monitor and predict fatigue crack propagation is presented in this work. The technique makes use of sequential Monte Carlo sampling combined with the probability density functions of the model parameters. The technique leads to an adaptive dynamic state-space model for crack evolution able to identify the most probable parameters, enhancing the estimation of the residual life of the system. The lifetime predictor presented here could be implemented in advanced maintenance strategies for critical structures employed in civil, industrial, and aerospace fields. The algorithm is first applied to a simulated crack growth, and then to some experimental crack propagations from laboratory tests on helicopter panels. The applicability within on-line continuous monitoring systems is discussed at the end of the paper.

Improving the efficiency of Monte Carlo Bayesian calibration of hydrologic models via model pre-emption

Bayesian inference via Markov Chain Monte Carlo (MCMC) sampling and sequential Monte Carlo (SMC) sampling are popular methods for uncertainty analysis in hydrological modelling. However, application of these methodologies can incur significant computational costs. This study investigated using model pre-emption for improving the computational efficiency of MCMC and SMC samplers in the context of hydrological modelling. The proposed pre-emption strategy facilitates early termination of low-likelihood simulations and results in reduction of unnecessary simulation time steps. The proposed approach is incorporated into two samplers and applied to the calibration of three rainfall–runoff models. Results show that overall pre-emption savings range from 5 to 21%. Furthermore, results indicate that pre-emption savings are greatest during the pre-convergence ‘burn-in’ period (i.e., between 8 and 39%) and decrease as the algorithms converge towards high likelihood regions of parameter space. The observed savings are achieved with absolutely no change in the posterior set of parameters.

A simulation approach for change-points on phylogenetic trees.

We observe n sequences at each of m sites and assume that they have evolved from an ancestral sequence that forms the root of a binary tree of known topology and branch lengths, but the sequence states at internal nodes are unknown. The topology of the tree and branch lengths are the same for all sites, but the parameters of the evolutionary model can vary over sites. We assume a piecewise constant model for these parameters, with an unknown number of change-points and hence a transdimensional parameter space over which we seek to perform Bayesian inference. We propose two novel ideas to deal with the computational challenges of such inference. Firstly, we approximate the model based on the time machine principle: the top nodes of the binary tree (near the root) are replaced by an approximation of the true distribution; as more nodes are removed from the top of the tree, the cost of computing the likelihood is reduced linearly in n. The approach introduces a bias, which we investigate empirically. Secondly, we develop a particle marginal Metropolis-Hastings (PMMH) algorithm, that employs a sequential Monte Carlo (SMC) sampler and can use the first idea. Our time-machine PMMH algorithm copes well with one of the bottle-necks of standard computational algorithms: the transdimensional nature of the posterior distribution. The algorithm is implemented on simulated and real data examples, and we empirically demonstrate its potential to outperform competing methods based on approximate Bayesian computation (ABC) techniques.

Efficient sequential Monte Carlo sampling for continuous monitoring of a radiation situation

The monitoring of a radiation situation around a nuclear power plant is a demanding task due to the high uncertainty of all involved variables and limited availability of measurements from a sparse monitoring network. Assessment of the situation requires experienced specialists who may be unavailable during critical times. Our goal is to provide an automated method of instant radiation situation assessment that does not underestimate its uncertainty. We propose a state-space model based on an atmospheric dispersion model, local correction of a numerical weather model, and a temporal model of the released activity. This state-space model is highly nonlinear and evaluation of the likelihood function requires extensive numerical calculations. The sequential Monte Carlo method is one of the few options for estimating the state recursively. Since the simple bootstrap approach yields an extremely computationally demanding algorithm, we investigate the use of existing techniques for the design of a more efficient proposal density. We propose combining the Laplace approximation and various population Monte Carlo methods. Data from an existing monitoring network were used to calibrate relevant parts of the model. Performance of the methods in a real radiation emergency situation is evaluated in a simulated experiment due to the lack of real data. The proposed tailor-made proposal was found to be much more computationally efficient than previously published methods. The adaptive Monte Carlo methods thus represent a compelling computational approach for the evaluation of probabilistic environmental models. The data used and a Python implementation of the methods are available as supplementary material online.

Sequential Monte Carlo samplers for semi-linear inverse problems and application to magnetoencephalography

We discuss the use of a recent class of sequential Monte Carlo methods for solving inverse problems characterized by a semi-linear structure, i.e. where the data depend linearly on a subset of variables and nonlinearly on the remaining ones. In this type of problems, under proper Gaussian assumptions one can marginalize the linear variables. This means that the Monte Carlo procedure needs only to be applied to the nonlinear variables, while the linear ones can be treated analytically; as a result, the Monte Carlo variance and/or the computational cost decrease. We use this approach to solve the inverse problem of magnetoencephalography, with a multi-dipole model for the sources. Here, data depend nonlinearly on the number of sources and their locations, and depend linearly on their current vectors. The semi-analytic approach enables us to estimate the number of dipoles and their location from a whole time-series, rather than a single time point, while keeping a low computational cost.