New sequential Monte Carlo methods for nonlinear dynamic systems

In this paper we present several new sequential Monte Carlo (SMC) algorithms for online estimation (filtering) of nonlinear dynamic systems. SMC has been shown to be a powerful tool for dealing with complex dynamic systems. It sequentially generates Monte Carlo samples from a proposal distribution, adjusted by a set of importance weight with respect to a target distribution, to facilitate statistical inferences on the characteristic (state) of the system. The key to a successful implementation of SMC in complex problems is the design of an efficient proposal distribution from which the Monte Carlo samples are generated. We propose several such proposal distributions that are efficient yet easy to generate samples from. They are efficient because they tend to utilize both the information in the state process and the observations. They are all Gaussian distributions hence are easy to sample from. The central ideas of the conventional nonlinear filters, such as extended Kalman filter, unscented Kalman filter and the Gaussian quadrature filter, are used to construct these proposal distributions. The effectiveness of the proposed algorithms are demonstrated through two applications—real time target tracking and the multiuser parameter tracking in CDMA communication systems.This work was supported in part by the U.S. National Science Foundation (NSF) under grants CCR-9875314, CCR-9980599, DMS-9982846, DMS-0073651 and DMS-0073601.     

Maximizing generalized linear mixed model likelihoods with an automated Monte Carlo EM algorithm

Two new implementations of the EM algorithm are proposed for maximum likelihood fitting of generalized linear mixed models. Both methods use random (independent and identically distributed) sampling to construct Monte Carlo approximations at the E-step. One approach involves generating random samples from the exact conditional distribution of the random effects (given the data) by rejection sampling, using the marginal distribution as a candidate. The second method uses a multivariate t importance sampling approximation. In many applications the two methods are complementary. Rejection sampling is more efficient when sample sizes are small, whereas importance sampling is better with larger sample sizes. Monte Carlo approximation using random samples allows the Monte Carlo error at each iteration to be assessed by using standard central limit theory combined with Taylor series methods. Specifically, we construct a sandwich variance estimate for the maximizer at each approximate E-step. This suggests a rule for automatically increasing the Monte Carlo sample size after iterations in which the true EM step is swamped by Monte Carlo error. In contrast, techniques for assessing Monte Carlo error have not been developed for use with alternative implementations of Monte Carlo EM algorithms utilizing Markov chain Monte Carlo E-step approximations. Three different data sets, including the infamous salamander data of McCullagh and Nelder, are used to illustrate the techniques and to compare them with the alternatives. The results show that the methods proposed can be considerably more efficient than those based on Markov chain Monte Carlo algorithms. However, the methods proposed may break down when the intractable integrals in the likelihood function are of high dimension.     

Approximate Methods for State-Space Models

State-space models provide an important body of techniques for analyzing time-series, but their use requires estimating unobserved states. The optimal estimate of the state is its conditional expectation given the observation histories, and computing this expectation is hard when there are nonlinearities. Existing filtering methods, including sequential Monte Carlo, tend to be either inaccurate or slow. In this paper, we study a nonlinear filter for nonlinear/non-Gaussian state-space models, which uses Laplace's method, an asymptotic series expansion, to approximate the state's conditional mean and variance, together with a Gaussian conditional distribution. This Laplace-Gaussian filter (LGF) gives fast, recursive, deterministic state estimates, with an error which is set by the stochastic characteristics of the model and is, we show, stable over time. We illustrate the estimation ability of the LGF by applying it to the problem of neural decoding and compare it to sequential Monte Carlo both in simulations and with real data. We find that the LGF can deliver superior results in a small fraction of the computing time.     

Nonparametric particle filtering and smoothing with quasi-Monte Carlo sampling

Sequential Monte Carlo methods (also known as particle filters and smoothers) are used for filtering and smoothing in general state-space models. These methods are based on importance sampling. In practice, it is often difficult to find a suitable proposal which allows effective importance sampling. This article develops an original particle filter and an original particle smoother which employ nonparametric importance sampling. The basic idea is to use a nonparametric estimate of the marginally optimal proposal. The proposed algorithms provide a better approximation of the filtering and smoothing distributions than standard methods. The methods’ advantage is most distinct in severely nonlinear situations. In contrast to most existing methods, they allow the use of quasi-Monte Carlo (QMC) sampling. In addition, they do not suffer from weight degeneration rendering a resampling step unnecessary. For the estimation of model parameters, an efficient on-line maximum-likelihood (ML) estimation technique is proposed which is also based on nonparametric approximations. All suggested algorithms have almost linear complexity for low-dimensional state-spaces. This is an advantage over standard smoothing and ML procedures. Particularly, all existing sequential Monte Carlo methods that incorporate QMC sampling have quadratic complexity. As an application, stochastic volatility estimation for high-frequency financial data is considered, which is of great importance in practice. The computer code is partly available as supplemental material.     

A pseudo-marginal sequential Monte Carlo algorithm for random effects models in Bayesian sequential design

Motivated by the need to sequentially design experiments for the collection of data in batches or blocks, a new pseudo-marginal sequential Monte Carlo algorithm is proposed for random effects models where the likelihood is not analytic, and has to be approximated. This new algorithm is an extension of the idealised sequential Monte Carlo algorithm where we propose to unbiasedly approximate the likelihood to yield an efficient exact-approximate algorithm to perform inference and make decisions within Bayesian sequential design. We propose four approaches to unbiasedly approximate the likelihood: standard Monte Carlo integration; randomised quasi-Monte Carlo integration, Laplace importance sampling and a combination of Laplace importance sampling and randomised quasi-Monte Carlo. These four methods are compared in terms of the estimates of likelihood weights and in the selection of the optimal sequential designs in an important pharmacological study related to the treatment of critically ill patients. As the approaches considered to approximate the likelihood can be computationally expensive, we exploit parallel computational architectures to ensure designs are derived in a timely manner.     

Some contributions to sequential Monte Carlo methods for option pricing

Pricing options is an important problem in financial engineering. In many scenarios of practical interest, financial option prices associated with an underlying asset reduces to computing an expectation w.r.t. a diffusion process. In general, these expectations cannot be calculated analytically, and one way to approximate these quantities is via the Monte Carlo (MC) method; MC methods have been used to price options since at least the 1970s. It has been seen in Del Moral P, Shevchenko PV. [Valuation of barrier options using sequential Monte Carlo. 2014. arXiv preprint] and Jasra A, Del Moral P. [Sequential Monte Carlo methods for option pricing. Stoch Anal Appl. 2011;29:292–316] that Sequential Monte Carlo (SMC) methods are a natural tool to apply in this context and can vastly improve over standard MC. In this article, in a similar spirit to Del Moral and Shevchenko (2014) and Jasra and Del Moral (2011), we show that one can achieve significant gains by using SMC methods by constructing a sequence of artificial target densities over time. In particular, we approximate the optimal importance sampling distribution in the SMC algorithm by using a sequence of weighting functions. This is demonstrated on two examples, barrier options and target accrual redemption notes (TARNs). We also provide a proof of unbiasedness of our SMC estimate.     

On spatiotemporal prediction for on-line monitoring data

Spatiotemporal prediction is of interest in many areas of applied statistics, especially in environmental monitoring with on-line data information. At first, this article reviews the approaches for spatiotemporal modeling in the context of stochastic processes and then introduces the new class of spatiotemporal dynamic linear models. Further, the methods for linear spatial data analysis, universal kriging and trend surface prediction, are related to the method of spatial linear Bayesian analysis. The Kalman filter is the preferred method for temporal linear Bayesian inferences. By combining the Kalman filter recursions with the trend surface predictor and universal kriging predictor, the prior and posterior spatiotemporal predictors for the observational process are derived, which form the main result of this article. The problem of spatiotemporal linear prediction in the case of unknown first and second order moments is treated as well.     

Polynomial Spline Estimation for a Generalized Additive Coefficient Model

Abstract.  We study a semiparametric generalized additive coefficient model (GACM), in which linear predictors in the conventional generalized linear models are generalized to unknown functions depending on certain covariates, and approximate the non-parametric functions by using polynomial spline. The asymptotic expansion with optimal rates of convergence for the estimators of the non-parametric part is established. Semiparametric generalized likelihood ratio test is also proposed to check if a non-parametric coefficient can be simplified as a parametric one. A conditional bootstrap version is suggested to approximate the distribution of the test under the null hypothesis. Extensive Monte Carlo simulation studies are conducted to examine the finite sample performance of the proposed methods. We further apply the proposed model and methods to a data set from a human visceral Leishmaniasis study conducted in Brazil from 1994 to 1997. Numerical results outperform the traditional generalized linear model and the proposed GACM is preferable.     

Hamiltonian Monte Carlo acceleration using surrogate functions with random bases

For big data analysis, high computational cost for Bayesian methods often limits their applications in practice. In recent years, there have been many attempts to improve computational efficiency of Bayesian inference. Here we propose an efficient and scalable computational technique for a state-of-the-art Markov chain Monte Carlo methods, namely, Hamiltonian Monte Carlo. The key idea is to explore and exploit the structure and regularity in parameter space for the underlying probabilistic model to construct an effective approximation of its geometric properties. To this end, we build a surrogate function to approximate the target distribution using properly chosen random bases and an efficient optimization process. The resulting method provides a flexible, scalable, and efficient sampling algorithm, which converges to the correct target distribution. We show that by choosing the basis functions and optimization process differently, our method can be related to other approaches for the construction of surrogate functions such as generalized additive models or Gaussian process models. Experiments based on simulated and real data show that our approach leads to substantially more efficient sampling algorithms compared to existing state-of-the-art methods.     

Fitting finite mixture models using iterative Monte Carlo classification

Parameters of a finite mixture model are often estimated by the expectation–maximization (EM) algorithm where the observed data log-likelihood function is maximized. This paper proposes an alternative approach for fitting finite mixture models. Our method, called the iterative Monte Carlo classification (IMCC), is also an iterative fitting procedure. Within each iteration, it first estimates the membership probabilities for each data point, namely the conditional probability of a data point belonging to a particular mixing component given that the data point value is obtained, it then classifies each data point into a component distribution using the estimated conditional probabilities and the Monte Carlo method. It finally updates the parameters of each component distribution based on the classified data. Simulation studies were conducted to compare IMCC with some other algorithms for fitting mixture normal, and mixture t, densities.     

Stopping-time resampling for sequential Monte Carlo methods

Summary.  Motivated by the statistical inference problem in population genetics, we present a new sequential importance sampling with resampling strategy. The idea of resampling is key to the recent surge of popularity of sequential Monte Carlo methods in the statistics and engin-eering communities, but existing resampling techniques do not work well for coalescent-based inference problems in population genetics. We develop a new method called 'stopping-time resampling', which allows us to compare partially simulated samples at different stages to terminate unpromising partial samples and to multiply promising samples early on. To illustrate the idea, we first apply the new method to approximate the solution of a Dirichlet problem and the likelihood function of a non-Markovian process. Then we focus on its application in population genetics. All our examples show that the new resampling method can significantly improve the computational efficiency of existing sequential importance sampling methods.     

Uncertainty Quantification in the Ensemble Kalman Filter

The ensemble Kalman filter (EnKF) provides an approximate, sequential Monte Carlo solution to the recursive data assimilation algorithm for hidden Markov chains. The challenging conditioning step is approximated by a linear updating, and the updating weights, termed Kalman weights, are inferred from the ensemble members. The EnKF scheme is known to provide unstable predictions and to underestimate the prediction intervals, and even sometimes to diverge. The underlying cause for these shortcomings is poorly understood. We find that the ensemble members couple in the conditioning procedure and that the coupling increase multiplicatively in the recursive conditioning steps. Under reasonable Gauss‐independence assumptions, exact expressions for this correlation are developed. Moreover, expressions for the precision of the predictions and the downward bias in the empirical variance introduced in one conditioning step are found. These results are confirmed by a Gauss‐linear simulation study. Furthermore, we quantitatively evaluate an alternative, improved EnKF scheme on the basis of transformations of ensemble members under the same Gauss‐independent assumptions. The scheme is compared with the frequently used ensemble inflation scheme.     

Importance sampling type estimators based on approximate marginal Markov chain Monte Carlo

We consider importance sampling (IS) type weighted estimators based on Markov chain Monte Carlo (MCMC) targeting an approximate marginal of the target distribution. In the context of Bayesian latent variable models, the MCMC typically operates on the hyperparameters, and the subsequent weighting may be based on IS or sequential Monte Carlo (SMC), but allows for multilevel techniques as well. The IS approach provides a natural alternative to delayed acceptance (DA) pseudo-marginal/particle MCMC, and has many advantages over DA, including a straightforward parallelization and additional flexibility in MCMC implementation. We detail minimal conditions which ensure strong consistency of the suggested estimators, and provide central limit theorems with expressions for asymptotic variances. We demonstrate how our method can make use of SMC in the state space models context, using Laplace approximations and time-discretized diffusions. Our experimental results are promising and show that the IS-type approach can provide substantial gains relative to an analogous DA scheme, and is often competitive even without parallelization.     

The Gibbs sampler with particle efficient importance sampling for state-space models*

We consider Particle Gibbs (PG) for Bayesian analysis of non-linear non-Gaussian state-space models. As a Monte Carlo (MC) approximation of the Gibbs procedure, PG uses sequential MC (SMC) importance sampling inside the Gibbs to update the latent states. We propose to combine PG with the Particle Efficient Importance Sampling (PEIS). By using SMC sampling densities which are approximately globally fully adapted to the targeted density of the states, PEIS can substantially improve the simulation efficiency of the PG relative to existing PG implementations. The efficiency gains are illustrated in PG applications to a non-linear local-level model and stochastic volatility models.     

On sequential Monte Carlo, partial rejection control and approximate Bayesian computation

We present a variant of the sequential Monte Carlo sampler by incorporating the partial rejection control mechanism of Liu (2001). We show that the resulting algorithm can be considered as a sequential Monte Carlo sampler with a modified mutation kernel. We prove that the new sampler can reduce the variance of the incremental importance weights when compared with standard sequential Monte Carlo samplers, and provide a central limit theorem. Finally, the sampler is adapted for application under the challenging approximate Bayesian computation modelling framework.     

On-line inference for hidden Markov models via particle filters

Summary.  We consider the on-line Bayesian analysis of data by using a hidden Markov model, where inference is tractable conditional on the history of the state of the hidden component. A new particle filter algorithm is introduced and shown to produce promising results when analysing data of this type. The algorithm is similar to the mixture Kalman filter but uses a different resampling algorithm. We prove that this resampling algorithm is computationally efficient and optimal, among unbiased resampling algorithms, in terms of minimizing a squared error loss function. In a practical example, that of estimating break points from well-log data, our new particle filter outperforms two other particle filters, one of which is the mixture Kalman filter, by between one and two orders of magnitude.     

On population-based simulation for static inference

In this paper we present a review of population-based simulation for static inference problems. Such methods can be described as generating a collection of random variables {X _n } _n=1,…,N in parallel in order to simulate from some target density π (or potentially sequence of target densities). Population-based simulation is important as many challenging sampling problems in applied statistics cannot be dealt with successfully by conventional Markov chain Monte Carlo (MCMC) methods. We summarize population-based MCMC (Geyer, Computing Science and Statistics: The 23rd Symposium on the Interface, pp. 156–163, 1991; Liang and Wong, J. Am. Stat. Assoc. 96, 653–666, 2001) and sequential Monte Carlo samplers (SMC) (Del Moral, Doucet and Jasra, J. Roy. Stat. Soc. Ser. B 68, 411–436, 2006a), providing a comparison of the approaches. We give numerical examples from Bayesian mixture modelling (Richardson and Green, J. Roy. Stat. Soc. Ser. B 59, 731–792, 1997).     

A SEMIPARAMETRIC MIXTURE MODEL FOR THE ANALYSIS OF COMPETING RISKS DATA

This paper deals with the regression analysis of failure time data when there are censoring and multiple types of failures. We propose a semiparametric generalization of a parametric mixture model of Larson & Dinse (1985), for which the marginal probabilities of the various failure types are logistic functions of the covariates. Given the type of failure, the conditional distribution of the time to failure follows a proportional hazards model. A marginal like lihood approach to estimating regression parameters is suggested, whereby the baseline hazard functions are eliminated as nuisance parameters. The Monte Carlo method is used to approximate the marginal likelihood; the resulting function is maximized easily using existing software. Some guidelines for choosing the number of Monte Carlo replications are given. Fixing the regression parameters at their estimated values, the full likelihood is maximized via an EM algorithm to estimate the baseline survivor functions. The methods suggested are illustrated using the Stanford heart transplant data.     

Inequality Constrained State-Space Models

ABSTRACT

The standard Kalman filter cannot handle inequality constraints imposed on the state variables, as state truncation induces a nonlinear and non-Gaussian model. We propose a Rao-Blackwellized particle filter with the optimal importance function for forward filtering and the likelihood function evaluation. The particle filter effectively enforces the state constraints when the Kalman filter violates them. Monte Carlo experiments demonstrate excellent performance of the proposed particle filter with Rao-Blackwellization, in which the Gaussian linear sub-structure is exploited at both the cross-sectional and temporal levels.     

Layered adaptive importance sampling

Monte Carlo methods represent the de facto standard for approximating complicated integrals involving multidimensional target distributions. In order to generate random realizations from the target distribution, Monte Carlo techniques use simpler proposal probability densities to draw candidate samples. The performance of any such method is strictly related to the specification of the proposal distribution, such that unfortunate choices easily wreak havoc on the resulting estimators. In this work, we introduce a layered (i.e., hierarchical) procedure to generate samples employed within a Monte Carlo scheme. This approach ensures that an appropriate equivalent proposal density is always obtained automatically (thus eliminating the risk of a catastrophic performance), although at the expense of a moderate increase in the complexity. Furthermore, we provide a general unified importance sampling (IS) framework, where multiple proposal densities are employed and several IS schemes are introduced by applying the so-called deterministic mixture approach. Finally, given these schemes, we also propose a novel class of adaptive importance samplers using a population of proposals, where the adaptation is driven by independent parallel or interacting Markov chain Monte Carlo (MCMC) chains. The resulting algorithms efficiently combine the benefits of both IS and MCMC methods.