Source author record

Ajay Jasra

Ajay Jasra appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record

Catalog footprint

What is connected

62works
16topics
4close collaborators

Actions

Connect this record

Log in to claim

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

62 published item(s)

preprint2026arXiv

New Trends in the Stability of Sinkhorn Semigroups

Entropic optimal transport problems play an increasingly important role in machine learning and generative modelling. In contrast with optimal transport maps which often have limited applicability in high dimensions, Schrodinger bridges can be solved using the celebrated Sinkhorn's algorithm, a.k.a. the iterative proportional fitting procedure. The stability properties of Sinkhorn bridges when the number of iterations tends to infinity is a very active research area in applied probability and machine learning. Traditional proofs of convergence are mainly based on nonlinear versions of Perron-Frobenius theory and related Hilbert projective metric techniques, gradient descent, Bregman divergence techniques and Hamilton-Jacobi-Bellman equations, including propagation of convexity profiles based on coupling diffusions by reflection methods. The objective of this review article is to present, in a self-contained manner, recently developed Sinkhorn/Gibbs-type semigroup analysis based upon contraction coefficients and Lyapunov-type operator-theoretic techniques. These powerful, off-the-shelf semigroup methods are based upon transportation cost inequalities (e.g. log-Sobolev, Talagrand quadratic inequality, curvature estimates), $ϕ$-divergences, Kantorovich-type criteria and Dobrushin contraction-type coefficients on weighted Banach spaces as well as Wasserstein distances. This novel semigroup analysis allows one to unify and simplify many arguments in the stability of Sinkhorn algorithm. It also yields new contraction estimates w.r.t. generalized $ϕ$-entropies, as well as weighted total variation norms, Kantorovich criteria and Wasserstein distances.

preprint2023arXiv

Antithetic Multilevel Particle Filters

In this paper we consider the filtering of partially observed multi-dimensional diffusion processes that are observed regularly at discrete times. This is a challenging problem which requires the use of advanced numerical schemes based upon time-discretization of the diffusion process and then the application of particle filters. Perhaps the state-of-the-art method for moderate dimensional problems is the multilevel particle filter of \cite{mlpf}. This is a method that combines multilevel Monte Carlo and particle filters. The approach in that article is based intrinsically upon an Euler discretization method. We develop a new particle filter based upon the antithetic truncated Milstein scheme of \cite{ml_anti}. We show that for a class of diffusion problems, for $ε>0$ given, that the cost to produce a mean square error (MSE) in estimation of the filter, of $\mathcal{O}(ε^2)$ is $\mathcal{O}(ε^{-2}\log(ε)^2)$. In the case of multidimensional diffusions with non-constant diffusion coefficient, the method of \cite{mlpf} has a cost of $\mathcal{O}(ε^{-2.5})$ to achieve the same MSE. We support our theory with numerical results in several examples.

preprint2022arXiv

A Lagged Particle Filter for Stable Filtering of certain High-Dimensional State-Space Models

We consider the problem of high-dimensional filtering of state-space models (SSMs) at discrete times. This problem is particularly challenging as analytical solutions are typically not available and many numerical approximation methods can have a cost that scales exponentially with the dimension of the hidden state. Inspired by lag-approximation methods for the smoothing problem, we introduce a lagged approximation of the smoothing distribution that is necessarily biased. For certain classes of SSMs, particularly those that forget the initial condition exponentially fast in time, the bias of our approximation is shown to be uniformly controlled in the dimension and exponentially small in time. We develop a sequential Monte Carlo (SMC) method to recursively estimate expectations with respect to our biased filtering distributions. Moreover, we prove for a class of class of SSMs that can contain dependencies amongst coordinates that as the dimension $d\rightarrow\infty$ the cost to achieve a stable mean square error in estimation, for classes of expectations, is of $\mathcal{O}(Nd^2)$ per-unit time, where $N$ is the number of simulated samples in the SMC algorithm. Our methodology is implemented on several challenging high-dimensional examples including the conservative shallow-water model.

preprint2022arXiv

Convergence Speed and Approximation Accuracy of Numerical MCMC

When implementing Markov Chain Monte Carlo (MCMC) algorithms, perturbation caused by numerical errors is sometimes inevitable. This paper studies how perturbation of MCMC affects the convergence speed and Monte Carlo estimation accuracy. Our results show that when the original Markov chain converges to stationarity fast enough and the perturbed transition kernel is a good approximation to the original transition kernel, the corresponding perturbed sampler has similar convergence speed and high approximation accuracy as well. We discuss two different analysis frameworks: ergodicity and spectral gap, both are widely used in the literature. Our results can be easily extended to obtain non-asymptotic error bounds for MCMC estimators. We also demonstrate how to apply our convergence and approximation results to the analysis of specific sampling algorithms, including Random walk Metropolis and Metropolis adjusted Langevin algorithm with perturbed target densities, and parallel tempering Monte Carlo with perturbed densities. Finally we present some simple numerical examples to verify our theoretical claims.

preprint2022arXiv

Unbiased Estimation of the Vanilla and Deterministic Ensemble Kalman-Bucy Filters

In this article we consider the development of an unbiased estimator for the ensemble Kalman--Bucy filter (EnKBF). The EnKBF is a continuous-time filtering methodology which can be viewed as a continuous-time analogue of the famous discrete-time ensemble Kalman filter. Our unbiased estimators will be motivated from recent work [Rhee \& Glynn 2010, [31]] which introduces randomization as a means to produce unbiased and finite variance estimators. The randomization enters through both the level of discretization, and through the number of samples at each level. Our estimator will be specific to linear and Gaussian settings, where we know that the EnKBF is consistent, in the particle limit $N \rightarrow \infty$, with the KBF. We highlight this for two particular variants of the EnKBF, i.e. the deterministic and vanilla variants, and demonstrate this on a linear Ornstein--Uhlenbeck process. We compare this with the EnKBF and the multilevel (MLEnKBF), for experiments with varying dimension size. We also provide a proof of the multilevel deterministic EnKBF, which provides a guideline for some of the unbiased methods.

preprint2022arXiv

Unbiased Parameter Inference for a Class of Partially Observed Levy-Process Models

We consider the problem of static Bayesian inference for partially observed Levy-process models. We develop a methodology which allows one to infer static parameters and some states of the process, without a bias from the time-discretization of the afore-mentioned Levy process. The unbiased method is exceptionally amenable to parallel implementation and can be computationally efficient relative to competing approaches. We implement the method on S & P 500 log-return daily data and compare it to some Markov chain Monte Carlo (MCMC) algorithm.

preprint2021arXiv

A 4D-Var Method with Flow-Dependent Background Covariances for the Shallow-Water Equations

The 4D-Var method for filtering partially observed nonlinear chaotic dynamical systems consists of finding the maximum a-posteriori (MAP) estimator of the initial condition of the system given observations over a time window, and propagating it forward to the current time via the model dynamics. This method forms the basis of most currently operational weather forecasting systems. In practice the optimization becomes infeasible if the time window is too long due to the non-convexity of the cost function, the effect of model errors, and the limited precision of the ODE solvers. Hence the window has to be kept sufficiently short, and the observations in the previous windows can be taken into account via a Gaussian background (prior) distribution. The choice of the background covariance matrix is an important question that has received much attention in the literature. In this paper, we define the background covariances in a principled manner, based on observations in the previous $b$ assimilation windows, for a parameter $b\ge 1$. The method is at most $b$ times more computationally expensive than using fixed background covariances, requires little tuning, and greatly improves the accuracy of 4D-Var. As a concrete example, we focus on the shallow-water equations. The proposed method is compared against state-of-the-art approaches in data assimilation and is shown to perform favourably on simulated data. We also illustrate our approach on data from the recent tsunami of 2011 in Fukushima, Japan.

preprint2021arXiv

Log-Normalization Constant Estimation using the Ensemble Kalman-Bucy Filter with Application to High-Dimensional Models

In this article we consider the estimation of the log-normalization constant associated to a class of continuous-time filtering models. In particular, we consider ensemble Kalman-Bucy filter based estimates based upon several nonlinear Kalman-Bucy diffusions. Based upon new conditional bias results for the mean of the afore-mentioned methods, we analyze the empirical log-scale normalization constants in terms of their $\mathbb{L}_n-$errors and conditional bias. Depending on the type of nonlinear Kalman-Bucy diffusion, we show that these are of order $(\sqrt{t/N}) + t/N$ or $1/\sqrt{N}$ ($\mathbb{L}_n-$errors) and of order $[t+\sqrt{t}]/N$ or $1/N$ (conditional bias), where $t$ is the time horizon and $N$ is the ensemble size. Finally, we use these results for online static parameter estimation for above filtering models and implement the methodology for both linear and nonlinear models.

preprint2021arXiv

On Unbiased Estimation for Discretized Models

In this article, we consider computing expectations w.r.t. probability measures which are subject to discretization error. Examples include partially observed diffusion processes or inverse problems, where one may have to discretize time and/or space, in order to practically work with the probability of interest. Given access only to these discretizations, we consider the construction of unbiased Monte Carlo estimators of expectations w.r.t. such target probability distributions. It is shown how to obtain such estimators using a novel adaptation of randomization schemes and Markov simulation methods. Under appropriate assumptions, these estimators possess finite variance and finite expected cost. There are two important consequences of this approach: (i) unbiased inference is achieved at the canonical complexity rate, and (ii) the resulting estimators can be generated independently, thereby allowing strong scaling to arbitrarily many parallel processors. Several algorithms are presented, and applied to some examples of Bayesian inference problems, with both simulated and real observed data.

preprint2021arXiv

Unbiased inference for discretely observed hidden Markov model diffusions

We develop a Bayesian inference method for diffusions observed discretely and with noise, which is free of discretisation bias. Unlike existing unbiased inference methods, our method does not rely on exact simulation techniques. Instead, our method uses standard time-discretised approximations of diffusions, such as the Euler--Maruyama scheme. Our approach is based on particle marginal Metropolis--Hastings, a particle filter, randomised multilevel Monte Carlo, and importance sampling type correction of approximate Markov chain Monte Carlo. The resulting estimator leads to inference without a bias from the time-discretisation as the number of Markov chain iterations increases. We give convergence results and recommend allocations for algorithm inputs. Our method admits a straightforward parallelisation, and can be computationally efficient. The user-friendly approach is illustrated on three examples, where the underlying diffusion is an Ornstein--Uhlenbeck process, a geometric Brownian motion, and a 2d non-reversible Langevin equation.

preprint2020arXiv

A practical and efficient approach for Bayesian quantum state estimation

Bayesian inference is a powerful paradigm for quantum state tomography, treating uncertainty in meaningful and informative ways. Yet the numerical challenges associated with sampling from complex probability distributions hampers Bayesian tomography in practical settings. In this Article, we introduce an improved, self-contained approach for Bayesian quantum state estimation. Leveraging advances in machine learning and statistics, our formulation relies on highly efficient preconditioned Crank--Nicolson sampling and a pseudo-likelihood. We theoretically analyze the computational cost, and provide explicit examples of inference for both actual and simulated datasets, illustrating improved performance with respect to existing approaches.

preprint2020arXiv

A Wasserstein Coupled Particle Filter for Multilevel Estimation

In this paper, we consider the filtering problem for partially observed diffusions, which are regularly observed at discrete times. We are concerned with the case when one must resort to time-discretization of the diffusion process if the transition density is not available in an appropriate form. In such cases, one must resort to advanced numerical algorithms such as particle filters to consistently estimate the filter. It is also well known that the particle filter can be enhanced by considering hierarchies of discretizations and the multilevel Monte Carlo (MLMC) method, in the sense of reducing the computational effort to achieve a given mean square error (MSE). A variety of multilevel particle filters (MLPF) have been suggested in the literature, e.g., in Jasra et al., SIAM J, Numer. Anal., 55, 3068--3096. Here we introduce a new alternative that involves a resampling step based on the optimal Wasserstein coupling. We prove a central limit theorem (CLT) for the new method. On considering the asymptotic variance, we establish that in some scenarios, there is a reduction, relative to the approach in the aforementioned paper by Jasra et al., in computational effort to achieve a given MSE. These findings are confirmed in numerical examples. We also consider filtering diffusions with unstable dynamics; we empirically show that in such cases a change of measure technique seems to be required to maintain our findings.

preprint2020arXiv

Multi-Index Sequential Monte Carlo Methods for partially observed Stochastic Partial Differential Equations

In this paper we consider sequential joint state and static parameter estimation given discrete time observations associated to a partially observed stochastic partial differential equation (SPDE). It is assumed that one can only estimate the hidden state using a discretization of the model. In this context, it is known that the multi-index Monte Carlo (MIMC) method of [11] can be used to improve over direct Monte Carlo from the most precise discretizaton. However, in the context of interest, it cannot be directly applied, but rather must be used within another advanced method such as sequential Monte Carlo (SMC). We show how one can use the MIMC method by renormalizing the MI identity and approximating the resulting identity using the SMC$^2$ method of [5]. We prove that our approach can reduce the cost to obtain a given mean square error (MSE), relative to just using SMC$^2$ on the most precise discretization. We demonstrate this with some numerical examples.

preprint2020arXiv

Multilevel Particle Filters for the Non-Linear Filtering Problem in Continuous Time

In the following article we consider the numerical approximation of the non-linear filter in continuous-time, where the observations and signal follow diffusion processes. Given access to high-frequency, but discrete-time observations, we resort to a first order time discretization of the non-linear filter, followed by an Euler discretization of the signal dynamics. In order to approximate the associated discretized non-linear filter, one can use a particle filter (PF). Under assumptions, this can achieve a mean square error of $\mathcal{O}(ε^2)$, for $ε>0$ arbitrary, such that the associated cost is $\mathcal{O}(ε^{-4})$. We prove, under assumptions, that the multilevel particle filter (MLPF) of Jasra et al (2017) can achieve a mean square error of $\mathcal{O}(ε^2)$, for cost $\mathcal{O}(ε^{-3})$. This is supported by numerical simulations in several examples.

preprint2020arXiv

Unbiased Estimation of the Gradient of the Log-Likelihood in Inverse Problems

We consider the problem of estimating a parameter associated to a Bayesian inverse problem. Treating the unknown initial condition as a nuisance parameter, typically one must resort to a numerical approximation of gradient of the log-likelihood and also adopt a discretization of the problem in space and/or time. We develop a new methodology to unbiasedly estimate the gradient of the log-likelihood with respect to the unknown parameter, i.e. the expectation of the estimate has no discretization bias. Such a property is not only useful for estimation in terms of the original stochastic model of interest, but can be used in stochastic gradient algorithms which benefit from unbiased estimates. Under appropriate assumptions, we prove that our estimator is not only unbiased but of finite variance. In addition, when implemented on a single processor, we show that the cost to achieve a given level of error is comparable to multilevel Monte Carlo methods, both practically and theoretically. However, the new algorithm provides the possibility for parallel computation on arbitrarily many processors without any loss of efficiency, asymptotically. In practice, this means any precision can be achieved in a fixed, finite constant time, provided that enough processors are available.

preprint2020arXiv

Unbiased Estimation of the Solution to Zakai's Equation

In the following article we consider the non-linear filtering problem in continuous-time and in particular the solution to Zakai's equation or the normalizing constant. We develop a methodology to produce finite variance, almost surely unbiased estimators of the solution to Zakai's equation. That is, given access to only a first order discretization of solution to the Zakai equation, we present a method which can remove this discretization bias. The approach, under assumptions, is proved to have finite variance and is numerically compared to using a particular multilevel Monte Carlo method.

preprint2020arXiv

Unbiased Filtering of a Class of Partially Observed Diffusions

In this article we consider a Monte Carlo-based method to filter partially observed diffusions observed at regular and discrete times. Given access only to Euler discretizations of the diffusion process, we present a new procedure which can return online estimates of the filtering distribution with no discretization bias and finite variance. Our approach is based upon a novel double application of the randomization methods of Rhee & Glynn (2015) along with the multilevel particle filter (MLPF) approach of Jasra et al (2017). A numerical comparison of our new approach with the MLPF, on a single processor, shows that similar errors are possible for a mild increase in computational cost. However, the new method scales strongly to arbitrarily many processors.

preprint2020arXiv

Uncertainty modelling and computational aspects of data association

A novel solution to the smoothing problem for multi-object dynamical systems is proposed and evaluated. The systems of interest contain an unknown and varying number of dynamical objects that are partially observed under noisy and corrupted observations. An alternative representation of uncertainty is considered in order to account for the lack of information about the different aspects of this type of complex system. The corresponding statistical model can be formulated as a hierarchical model consisting of conditionally-independent hidden Markov models. This particular structure is leveraged to propose an efficient method in the context of Markov chain Monte Carlo (MCMC) by relying on an approximate solution to the corresponding filtering problem, in a similar fashion to particle MCMC. This approach is shown to outperform existing algorithms in a range of scenarios.

preprint2018arXiv

Central Limit Theorems for Coupled Particle Filters

In this article we prove a new central limit theorem (CLT) for coupled particle filters (CPFs). CPFs are used for the sequential estimation of the difference of expectations w.r.t. filters which are in some sense close. Examples include the estimation of the filtering distribution associated to different parameters (finite difference estimation) and filters associated to partially observed discretized diffusion processes (PODDP) and the implementation of the multilevel Monte Carlo (MLMC) identity. We develop new theory for CPFs and based upon several results, we propose a new CPF which approximates the maximal coupling (MCPF) of a pair of predictor distributions. In the context of ML estimation associated to PODDP with discretization $Δ_l$ we show that the MCPF and the approach in Jasra et al. (2018) have, under assumptions, an asymptotic variance that is upper-bounded by an expression that is (almost) $\mathcal{O}(Δ_l)$, uniformly in time. The $\mathcal{O}(Δ_l)$ rate preserves the so-called forward rate of the diffusion in some scenarios which is not the case for the CPF in Jasra et al (2017).

preprint2016arXiv

A Note on Random Walks with Absorbing barriers and Sequential Monte Carlo Methods

In this article we consider importance sampling (IS) and sequential Monte Carlo (SMC) methods in the context of 1-dimensional random walks with absorbing barriers. In particular, we develop a very precise variance analysis for several IS and SMC procedures. We take advantage of some explicit spectral formulae available for these models to derive sharp and explicit estimates; this provides stability properties of the associated normalized Feynman-Kac semigroups. Our analysis allows one to compare the variance of SMC and IS techniques for these models. The work in this article, is one of the few to consider an in-depth analysis of an SMC method for a particular model-type as well as variance comparison of SMC algorithms.

preprint2016arXiv

Flexible online multivariate regression with variational Bayes and the matrix-variate Dirichlet process

Flexible regression methods where interest centres on the way that the whole distribution of a response vector changes with covariates are very useful in some applications. A recently developed technique in this regard uses the matrix-variate Dirichlet process as a prior for a mixing distribution on a coefficient in a multivariate linear regression model. The method is attractive, particularly in the multivariate setting, for the convenient way that it allows for borrowing strength across different component regressions and for its computational simplicity and tractability. The purpose of the present article is to develop fast online variational Bayes approaches to fitting this model and to investigate how they perform compared to MCMC and batch variational methods in a number of scenarios.

preprint2016arXiv

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.

preprint2016arXiv

Multilevel Particle Filters: Normalizing Constant Estimation

In this article we introduce two new estimates of the normalizing constant (or marginal likelihood) for partially observed diffusion (POD) processes, with discrete observations. One estimate is biased but non-negative and the other is unbiased but not almost surely non-negative. Our method uses the multilevel particle filter of Jasra et al (2015). We show that, under assumptions, for Euler discretized PODs and a given $\varepsilon>0$. in order to obtain a mean square error (MSE) of $\mathcal{O}(\varepsilon^2)$ one requires a work of $\mathcal{O}(\varepsilon^{-2.5})$ for our new estimates versus a standard particle filter that requires a work of $\mathcal{O}(\varepsilon^{-3})$. Our theoretical results are supported by numerical simulations.

preprint2016arXiv

Multilevel Sequential Monte Carlo Samplers for Normalizing Constants

This article considers the sequential Monte Carlo (SMC) approximation of ratios of normalizing constants associated to posterior distributions which in principle rely on continuum models. Therefore, the Monte Carlo estimation error and the discrete approximation error must be balanced. A multilevel strategy is utilized to substantially reduce the cost to obtain a given error level in the approximation as compared to standard estimators. Two estimators are considered and relative variance bounds are given. The theoretical results are numerically illustrated for the example of identifying a parametrized permeability in an elliptic equation given point-wise observations of the pressure.

preprint2016arXiv

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 to 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 method; Monte Carlo methods have been used to price options since at least the 1970's. It has been seen in Del Moral, P. \& Shevchenko, P.V. (2014) `Valuation of barrier options using Sequential Monte Carlo' and Jasra, A. \& Del Moral, P. (2011) `Sequential Monte Carlo for option pricing' that Sequential Monte Carlo (SMC) methods are a natural tool to apply in this context and can vastly improve over standard Monte Carlo. In this article, in a similar spirit to Del Moral, P. \& Shevchenko, P.V. (2014) `Valuation of barrier options using sequential Monte Carlo' and Jasra, A. \& Del Moral, P. (2011) `Sequential Monte Carlo for option pricing' 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 (TARN's). We also provide a proof of unbiasedness of our SMC estimate.

preprint2015arXiv

Bayesian Inference for Duplication-Mutation with Complementarity Network Models

We observe an undirected graph $G$ without multiple edges and self-loops, which is to represent a protein-protein interaction (PPI) network. We assume that $G$ evolved under the duplication-mutation with complementarity (DMC) model from a seed graph, $G_0$, and we also observe the binary forest $Γ$ that represents the duplication history of $G$. A posterior density for the DMC model parameters is established, and we outline a sampling strategy by which one can perform Bayesian inference; that sampling strategy employs a particle marginal Metropolis-Hastings (PMMH) algorithm. We test our methodology on numerical examples to demonstrate a high accuracy and precision in the inference of the DMC model's mutation and homodimerization parameters.

preprint2015arXiv

Biased Online Parameter Inference for State-Space Models

We consider Bayesian online static parameter estimation for state-space models. This is a very important problem, but is very computationally challenging as the state- of-the art methods that are exact, often have a computational cost that grows with the time parameter; perhaps the most successful algorithm is that of SMC2 [9]. We present a version of the SMC2 algorithm which has computational cost that does not grow with the time parameter. In addition, under assumptions, the algorithm is shown to provide consistent estimates of expectations w.r.t. the posterior. However, the cost to achieve this consistency can be exponential in the dimension of the parameter space; if this exponential cost is avoided, typically the algorithm is biased. The bias is investigated from a theoretical perspective and, under assumptions, we find that the bias does not accumulate as the time parameter grows. The algorithm is implemented on several Bayesian statistical models.

preprint2015arXiv

Multilevel particle filter

In this paper the filtering of partially observed diffusions, with discrete-time observations, is considered. It is assumed that only biased approximations of the diffusion can be obtained, for choice of an accuracy parameter indexed by $l$. A multilevel estimator is proposed, consisting of a telescopic sum of increment estimators associated to the successive levels. The work associated to $\mathcal{O}(\varepsilon^2)$ mean-square error between the multilevel estimator and average with respect to the filtering distribution is shown to scale optimally, for example as $\mathcal{O}(\varepsilon^{-2})$ for optimal rates of convergence of the underlying diffusion approximation. The method is illustrated on some toy examples as well as estimation of interest rate based on real S&P 500 stock price data.

preprint2014arXiv

A Sharp First Order Analysis of Feynman-Kac Particle Models

This article provides a new theory for the analysis of forward and backward particle approximations of Feynman-Kac models. Such formulae are found in a wide variety of applications and their numerical (particle) approximation are required due to their intractability. Under mild assumptions, we provide sharp and non-asymptotic first order expansions of these particle methods, potentially on path space and for possibly unbounded functions. These expansions allows one to consider upper and lower bound bias type estimates for a given time horizon $n$ and particle number $N$; these non-asymptotic estimates are of order $\mathcal{O}(n/N)$. Our approach is extended to tensor products of particle density profiles, leading to new sharp and non-asymptotic propagation of chaos estimates. The resulting upper and lower bound propagation of chaos estimates seems to be the first result of this kind for mean field particle models. As a by-product of our results, we also provide some analysis of the particle Gibbs sampler, providing first order expansions of the kernel and minorization estimates.

preprint2014arXiv

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 trans-dimensional 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 trans-dimensional 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.

preprint2014arXiv

A Stable Particle Filter in High-Dimensions

We consider the numerical approximation of the filtering problem in high dimensions, that is, when the hidden state lies in $\mathbb{R}^d$ with $d$ large. For low dimensional problems, one of the most popular numerical procedures for consistent inference is the class of approximations termed particle filters or sequential Monte Carlo methods. However, in high dimensions, standard particle filters (e.g. the bootstrap particle filter) can have a cost that is exponential in $d$ for the algorithm to be stable in an appropriate sense. We develop a new particle filter, called the \emph{space-time particle filter}, for a specific family of state-space models in discrete time. This new class of particle filters provide consistent Monte Carlo estimates for any fixed $d$, as do standard particle filters. Moreover, we expect that the state-space particle filter will scale much better with $d$ than the standard filter. We illustrate this analytically for a model of a simple i.i.d. structure and one of a Markovian structure in the $d$-dimensional space-direction, when we show that the algorithm exhibits certain stability properties as $d$ increases at a cost $\mathcal{O}(nNd^2)$, where $n$ is the time parameter and $N$ is the number of Monte Carlo samples, that are fixed and independent of $d$. Similar results are expected to hold, under a more general structure than the i.i.d.~one. independently of the dimension. Our theoretical results are also supported by numerical simulations on practical models of complex structures. The results suggest that it is indeed possible to tackle some high dimensional filtering problems using the space-time particle filter that standard particle filters cannot handle.

preprint2014arXiv

Approximate Bayesian Computation for a Class of Time Series Models

In the following article we consider approximate Bayesian computation (ABC) for certain classes of time series models. In particular, we focus upon scenarios where the likelihoods of the observations and parameter are intractable, by which we mean that one cannot evaluate the likelihood even up-to a positive unbiased estimate. This paper reviews and develops a class of approximation procedures based upon the idea of ABC, but, specifically maintains the probabilistic structure of the original statistical model. This idea is useful, in that it can facilitate an analysis of the bias of the approximation and the adaptation of established computational methods for parameter inference. Several existing results in the literature are surveyed and novel developments with regards to computation are given.

preprint2014arXiv

On the Convergence of Adaptive Sequential Monte Carlo Methods

In several implementations of Sequential Monte Carlo (SMC) methods it is natural, and important in terms of algorithmic efficiency, to exploit the information of the history of the samples to optimally tune their subsequent propagations. In this article we provide a carefully formulated asymptotic theory for a class of such \emph{adaptive} SMC methods. The theoretical framework developed here will cover, under assumptions, several commonly used SMC algorithms. There are only limited results about the theoretical underpinning of such adaptive methods: we will bridge this gap by providing a weak law of large numbers (WLLN) and a central limit theorem (CLT) for some of these algorithms. The latter seems to be the first result of its kind in the literature and provides a formal justification of algorithms used in many real data context. We establish that for a general class of adaptive SMC algorithms the asymptotic variance of the estimators from the adaptive SMC method is \emph{identical} to a so-called `perfect' SMC algorithm which uses ideal proposal kernels. Our results are supported by application on a complex high-dimensional posterior distribution associated with the Navier-Stokes model, where adapting high-dimensional parameters of the proposal kernels is critical for the efficiency of the algorithm.

preprint2014arXiv

Sequential Monte Carlo Methods for Bayesian Elliptic Inverse Problems

In this article we consider a Bayesian inverse problem associated to elliptic partial differential equations (PDEs) in two and three dimensions. This class of inverse problems is important in applications such as hydrology, but the complexity of the link function between unknown field and measurements can make it difficult to draw inference from the associated posterior. We prove that for this inverse problem a basic SMC method has a Monte Carlo rate of convergence with constants which are independent of the dimension of the discretization of the problem; indeed convergence of the SMC method is established in a function space setting. We also develop an enhancement of the sequential Monte Carlo (SMC) methods for inverse problems which were introduced in \cite{kantas}; the enhancement is designed to deal with the additional complexity of this elliptic inverse problem. The efficacy of the methodology, and its desirable theoretical properties, are demonstrated on numerical examples in both two and three dimensions.

preprint2014arXiv

Theory of Parallel Particle Filters for Hidden Markov Models

The objective of this article is to study the asymptotic behavior of a new particle filtering approach in the context of hidden Markov models (HMMs). In particular, we develop an algorithm where the latent-state sequence is segmented into multiple shorter portions, with an estimation technique based upon a separate particle filter in each portion. The partitioning facilitates the use of parallel processing. Based upon this approach, we introduce new estimators of the latent states and likelihood which have similar or better variance properties compared to estimators derived from standard particle filters. Moreover due to parallelization there is less wall-clock computational time. We show that the likelihood function estimator is unbiased, prove central limit theorem convergences of estimators, and provide consistent in-sample estimation of the asymptotic variances. The theoretical analyses, supported by a numerical study, show that the segmentation reduces the variances in smoothed latent-state estimators, in addition to the savings in wall-clock time.

preprint2013arXiv

A sequential algorithm for fast fitting of Dirichlet process mixture models

In this article we propose an improvement on the sequential updating and greedy search (SUGS) algorithm Wang and Dunson for fast fitting of Dirichlet process mixture models. The SUGS algorithm provides a means for very fast approximate Bayesian inference for mixture data which is particularly of use when data sets are so large that many standard Markov chain Monte Carlo (MCMC) algorithms cannot be applied efficiently, or take a prohibitively long time to converge. In particular, these ideas are used to initially interrogate the data, and to refine models such that one can potentially apply exact data analysis later on. SUGS relies upon sequentially allocating data to clusters and proceeding with an update of the posterior on the subsequent allocations and parameters which assumes this allocation is correct. Our modification softens this approach, by providing a probability distribution over allocations, with a similar computational cost; this approach has an interpretation as a variational Bayes procedure and hence we term it variational SUGS (VSUGS). It is shown in simulated examples that VSUGS can out-perform, in terms of density estimation and classification, the original SUGS algorithm in many scenarios. In addition, we present a data analysis for flow cytometry data, and SNP data via a three-class dirichlet process mixture model illustrating the apparent improvement over SUGS.

preprint2013arXiv

An Adaptive Sequential Monte Carlo Algorithm for Computing Permanents

We consider the computation of the permanent of a binary n by n matrix. It is well- known that the exact computation is a #P complete problem. A variety of Markov chain Monte Carlo (MCMC) computational algorithms have been introduced in the literature whose cost, in order to achieve a given level of accuracy, is O(n^7 log^4(n)). These algorithms use a particular collection of probability distributions, the `ideal' of which, (in some sense) are not known and need to be approximated. In this paper we propose an adaptive sequential Monte Carlo (SMC) algorithm that can both estimate the permanent and the ideal sequence of probabilities on the fly, with little user input. We provide theoretical results associated to the SMC estimate of the permanent, establishing its convergence and analyzing the relative variance of the estimate, in particular computating explicit bounds on the relative variance which depend upon n. Using this latter result, we provide a lower-bound on the computational cost, in order to achieve an arbitrarily small relative variance; we find that this cost is O(n^4 log^4(n)). Some numerical simulations are also given.

preprint2013arXiv

Approximate Inference for Observation Driven Time Series Models with Intractable Likelihoods

In the following article we consider approximate Bayesian parameter inference for observation driven time series models. Such statistical models appear in a wide variety of applications, including econometrics and applied mathematics. This article considers the scenario where the likelihood function cannot be evaluated point-wise; in such cases, one cannot perform exact statistical inference, including parameter estimation, which often requires advanced computational algorithms, such as Markov chain Monte Carlo (MCMC). We introduce a new approximation based upon approximate Bayesian computation (ABC). Under some conditions, we show that as $n\rightarrow\infty$, with $n$ the length of the time series, the ABC posterior has, almost surely, a maximum \emph{a posteriori} (MAP) estimator of the parameters which is different from the true parameter. However, a noisy ABC MAP, which perturbs the original data, asymptotically converges to the true parameter, almost surely. In order to draw statistical inference, for the ABC approximation adopted, standard MCMC algorithms can have acceptance probabilities that fall at an exponential rate in $n$ and slightly more advanced algorithms can mix poorly. We develop a new and improved MCMC kernel, which is based upon an exact approximation of a marginal algorithm, whose cost per-iteration is random but the expected cost, for good performance, is shown to be $\mathcal{O}(n^2)$ per-iteration.

preprint2013arXiv

Computational Methods for a Class of Network Models

In the following article we provide an exposition of exact computational methods to perform parameter inference from partially observed network models. In particular, we consider the duplication attachment (DA) model which has a likelihood function that typically cannot be evaluated in any reasonable computational time. We consider a number of importance sampling (IS) and sequential Monte Carlo (SMC) methods for approximating the likelihood of the network model for a fixed parameter value. It is well-known that for IS, the relative variance of the likelihood estimate typically grows at an exponential rate in the time parameter (here this is associated to the size of the network): we prove that, under assumptions, the SMC method will have relative variance which can grow only polynomially. In order to perform parameter estimation, we develop particle Markov chain Monte Carlo (PMCMC) algorithms to perform Bayesian inference. Such algorithms use the afore-mentioned SMC algorithms within the transition dynamics. The approaches are illustrated numerically.

preprint2013arXiv

On the Behaviour of the Backward Interpretation of Feynman-Kac Formulae under Verifiable Conditions

In the following article we consider the time-stability associated to the sequential Monte Carlo (SMC) estimate of the backward interpretation of Feynman-Kac Formulae. This is particularly of interest in the context of performing smoothing for hidden Markov models (HMMs). We prove a central limit theorem (CLT) under weaker assumptions than adopted in the literature. We then show that the associated asymptotic variance expression, for additive functionals grows at most linearly in time, under hypotheses that are weaker than those currently existing in the literature. The assumptions are verified for some state-space models.

preprint2013arXiv

Parameter Estimation in Hidden Markov Models with Intractable Likelihoods Using Sequential Monte Carlo

We propose sequential Monte Carlo based algorithms for maximum likelihood estimation of the static parameters in hidden Markov models with an intractable likelihood using ideas from approximate Bayesian computation. The static parameter estimation algorithms are gradient based and cover both offline and online estimation. We demonstrate their performance by estimating the parameters of three intractable models, namely the alpha-stable distribution, g-and-k distribution, and the stochastic volatility model with alpha-stable returns, using both real and synthetic data.

preprint2013arXiv

Sequential Monte Carlo Methods for High-Dimensional Inverse Problems: A case study for the Navier-Stokes equations

We consider the inverse problem of estimating the initial condition of a partial differential equation, which is only observed through noisy measurements at discrete time intervals. In particular, we focus on the case where Eulerian measurements are obtained from the time and space evolving vector field, whose evolution obeys the two-dimensional Navier-Stokes equations defined on a torus. This context is particularly relevant to the area of numerical weather forecasting and data assimilation. We will adopt a Bayesian formulation resulting from a particular regularization that ensures the problem is well posed. In the context of Monte Carlo based inference, it is a challenging task to obtain samples from the resulting high dimensional posterior on the initial condition. In real data assimilation applications it is common for computational methods to invoke the use of heuristics and Gaussian approximations. The resulting inferences are biased and not well-justified in the presence of non-linear dynamics and observations. On the other hand, Monte Carlo methods can be used to assimilate data in a principled manner, but are often perceived as inefficient in this context due to the high-dimensionality of the problem. In this work we will propose a generic Sequential Monte Carlo (SMC) sampling approach for high dimensional inverse problems that overcomes these difficulties. The method builds upon Markov chain Monte Carlo (MCMC) techniques, which are currently considered as benchmarks for evaluating data assimilation algorithms used in practice. In our numerical examples, the proposed SMC approach achieves the same accuracy as MCMC but in a much more efficient manner.

preprint2013arXiv

The Alive Particle Filter

In the following article we develop a particle filter for approximating Feynman-Kac models with indicator potentials. Examples of such models include approximate Bayesian computation (ABC) posteriors associated with hidden Markov models (HMMs) or rare-event problems. Such models require the use of advanced particle filter or Markov chain Monte Carlo (MCMC) algorithms e.g. Jasra et al. (2012), to perform estimation. One of the drawbacks of existing particle filters, is that they may 'collapse', in that the algorithm may terminate early, due to the indicator potentials. In this article, using a special case of the locally adaptive particle filter in Lee et al. (2013), which is closely related to Le Gland & Oudjane (2004), we use an algorithm which can deal with this latter problem, whilst introducing a random cost per-time step. This algorithm is investigated from a theoretical perspective and several results are given which help to validate the algorithms and to provide guidelines for their implementation. In addition, we show how this algorithm can be used within MCMC, using particle MCMC (Andrieu et al. 2010). Numerical examples are presented for ABC approximations of HMMs.

preprint2013arXiv

The TimeMachine for Inference on Stochastic Trees

The simulation of genealogical trees backwards in time, from observations up to the most recent common ancestor (MRCA), is hindered by the fact that, while approaching the root of the tree, coalescent events become rarer, with a corresponding increase in computation time. The recently proposed "Time Machine" tackles this issue by stopping the simulation of the tree before reaching the MRCA and correcting for the induced bias. We present a computationally efficient implementation of this approach that exploits multithreading.

preprint2013arXiv

Twisting the Alive Particle Filter

This work focuses on sampling from hidden Markov models (Cappe et al, 2005) whose observations have intractable density functions. We develop a new sequential Monte Carlo (Doucet et al, 2000 and Gordon et al, 1993) algorithm and a new particle marginal Metropolis-Hastings (Andrieu et al, 2010) algorithm for these purposes. We build from Jasra, et al (2013) and Whiteley, et al (2013) to construct the sequential Monte Carlo (SMC) algorithm (which we call the alive twisted particle filter). Like the alive particle filter of Jasra, et al (2013), our new SMC algorithm adopts an approximate Bayesian computation (Tavare et al, 1997) estimate of the HMM. Our alive twisted particle filter also uses a twisted proposal as in Whiteley, et al (2013) to obtain a low-variance estimate of the HMM normalising constant. We demonstrate via numerical examples that, in some scenarios, this estimate has a much lower variance than that of the estimate obtained via the alive particle filter. The low variance of this normalising constant estimate encourages the implementation of our SMC algorithm within a particle marginal Metropolis-Hastings (PMMH) scheme, and we call the resulting methodology ``alive twisted PMMH''. We numerically demonstrate on a stochastic volatility model how our alive twisted PMMH can converge faster than the standard alive PMMH of Jasra, et al (2013).

preprint2012arXiv

A Bayesian Mixture of Lasso Regressions with t-Errors

Motivated by a challenging problem in financial trading we are presented with a mixture of regressions with variable selection problem. In this regard, one is faced with data which possess outliers, skewness and, simultaneously, due to the nature of financial trading, one would like to be able to construct clusters with specific predictors that are fairly sparse. We develop a Bayesian mixture of lasso regressions with $t-$errors to reflect these specific demands. The resulting model is necessarily complex and to fit the model to real data, we develop a state-of-the-art Particle Markov chain Monte Carlo (PMCMC) algorithm based upon sequential Monte Carlo (SMC) methods. The model and algorithm are investigated on both simulated and real data.

preprint2012arXiv

Approximate Bayesian Computation for Smoothing

We consider a method for approximate inference in hidden Markov models (HMMs). The method circumvents the need to evaluate conditional densities of observations given the hidden states. It may be considered an instance of Approximate Bayesian Computation (ABC) and it involves the introduction of auxiliary variables valued in the same space as the observations. The quality of the approximation may be controlled to arbitrary precision through a parameter ε>0 . We provide theoretical results which quantify, in terms of ε, the ABC error in approximation of expectations of additive functionals with respect to the smoothing distributions. Under regularity assumptions, this error is O(nε), where n is the number of time steps over which smoothing is performed. For numerical implementation we adopt the forward-only sequential Monte Carlo (SMC) scheme of [16] and quantify the combined error from the ABC and SMC approximations. This forms some of the first quantitative results for ABC methods which jointly treat the ABC and simulation errors, with a finite number of data and simulated samples. When the HMM has unknown static parameters, we consider particle Markov chain Monte Carlo [2] (PMCMC) methods for batch statistical inference.

preprint2012arXiv

Bayesian Parameter Inference for Partially Observed Stopped Processes

In this article we consider Bayesian parameter inference associated to partially-observed stochastic processes that start from a set B0 and are stopped or killed at the first hitting time of a known set A. Such processes occur naturally within the context of a wide variety of applications. The associated posterior distributions are highly complex and posterior parameter inference requires the use of advanced Markov chain Monte Carlo (MCMC) techniques. Our approach uses a recently introduced simulation methodology, particle Markov chain Monte Carlo (PMCMC) (Andrieu et. al. 2010 [1]), where sequential Monte Carlo (SMC) approximations (see Doucet et. al. 2001 [18] and Liu 2001 [27]) are embedded within MCMC. However, when the parameter of interest is fixed, standard SMC algorithms are not always appropriate for many stopped processes. In Chen et. al. [11] and Del Moral 2004 [15] the authors introduce SMC approximations of multi-level Feynman-Kac formulae, which can lead to more efficient algorithms. This is achieved by devising a sequence of nested sets from B0 to A and then perform the resampling step only when the samples of the process reach intermediate level sets in the sequence. Naturally, the choice of the intermediate level sets is critical to the performance of such a scheme. In this paper, we demonstrate that multi-level SMC algorithms can be used as a proposal in PMCMC. In addition, we propose a flexible strategy that adapts the level sets for different parameter proposals. Our methodology is illustrated on the coalescent model with migration.

preprint2012arXiv

Inference for a Class of Partially Observed Point Process Models

This paper presents a simulation-based framework for sequential inference from partially and discretely observed point process (PP's) models with static parameters. Taking on a Bayesian perspective for the static parameters, we build upon sequential Monte Carlo (SMC) methods, investigating the problems of performing sequential filtering and smoothing in complex examples, where current methods often fail. We consider various approaches for approximating posterior distributions using SMC. Our approaches, with some theoretical discussion are illustrated on a doubly stochastic point process applied in the context of finance.

preprint2012arXiv

Linear Variance Bounds for Particle Approximations of Time-Homogeneous Feynman-Kac Formulae

This article establishes sufficient conditions for a linear-in-time bound on the non-asymptotic variance of particle approximations of time-homogeneous Feynman-Kac formulae. These formulae appear in a wide variety of applications including option pricing in finance and risk sensitive control in engineering. In direct Monte Carlo approximation of these formulae, the non-asymptotic variance typically increases at an exponential rate in the time parameter. It is shown that a linear bound holds when a non-negative kernel, defined by the logarithmic potential function and Markov kernel which specify the Feynman-Kac model, satisfies a type of multiplicative drift condition and other regularity assumptions. Examples illustrate that these conditions are general and flexible enough to accommodate two rather extreme cases, which can occur in the context of a non-compact state space: 1) when the potential function is bounded above, not bounded below and the Markov kernel is not ergodic; and 2) when the potential function is not bounded above, but the Markov kernel itself satisfies a multiplicative drift condition.

preprint2012arXiv

Marginal Likelihood Computation for Hidden Markov Models via Generalized Two-Filter Smoothing

In this note we introduce an estimate for the marginal likelihood associated to hidden Markov models (HMMs) using sequential Monte Carlo (SMC) approximations of the generalized two-filter smoothing decomposition (Briers, 2010). This estimate is shown to be unbiased and a central limit theorem (CLT) is established. This latter CLT also allows one to prove a CLT associated to estimates of expectations w.r.t. a marginal of the joint smoothing distribution; these form some of the first theoretical results associated to the SMC approximation of the generalized two-filter smoothing decomposition. The new estimate and its application is investigated from a numerical perspective.

preprint2012arXiv

On adaptive resampling strategies for sequential Monte Carlo methods

Sequential Monte Carlo (SMC) methods are a class of techniques to sample approximately from any sequence of probability distributions using a combination of importance sampling and resampling steps. This paper is concerned with the convergence analysis of a class of SMC methods where the times at which resampling occurs are computed online using criteria such as the effective sample size. This is a popular approach amongst practitioners but there are very few convergence results available for these methods. By combining semigroup techniques with an original coupling argument, we obtain functional central limit theorems and uniform exponential concentration estimates for these algorithms.

preprint2012arXiv

On the Stability of Sequential Monte Carlo Methods in High Dimensions

We investigate the stability of a Sequential Monte Carlo (SMC) method applied to the problem of sampling from a target distribution on $\mathbb{R}^d$ for large $d$. It is well known that using a single importance sampling step one produces an approximation for the target that deteriorates as the dimension $d$ increases, unless the number of Monte Carlo samples $N$ increases at an exponential rate in $d$. We show that this degeneracy can be avoided by introducing a sequence of artificial targets, starting from a `simple' density and moving to the one of interest, using an SMC method to sample from the sequence. Using this class of SMC methods with a fixed number of samples, one can produce an approximation for which the effective sample size (ESS) converges to a random variable $\varepsilon_N$ as $d\rightarrow\infty$ with $1<\varepsilon_{N}<N$. The convergence is achieved with a computational cost proportional to $Nd^2$. If $\varepsilon_N\ll N$, we can raise its value by introducing a number of resampling steps, say $m$ (where $m$ is independent of $d$). In this case, ESS converges to a random variable $\varepsilon_{N,m}$ as $d\rightarrow\infty$ and $\lim_{m\to\infty}\varepsilon_{N,m}=N$. Also, we show that the Monte Carlo error for estimating a fixed dimensional marginal expectation is of order $\frac{1}{\sqrt{N}}$ uniformly in $d$. The results imply that, in high dimensions, SMC algorithms can efficiently control the variability of the importance sampling weights and estimate fixed dimensional marginals at a cost which is less than exponential in $d$ and indicate that, in high dimensions, resampling leads to a reduction in the Monte Carlo error and increase in the ESS.

preprint2012arXiv

Robust model-based clustering with gene ranking

Cluster analysis of biological samples using gene expression measurements is a common task which aids the discovery of heterogeneous biological sub-populations having distinct mRNA profiles. Several model-based clustering algorithms have been proposed in which the distribution of gene expression values within each sub-group is assumed to be Gaussian. In the presence of noise and extreme observations, a mixture of Gaussian densities may over-fit and overestimate the true number of clusters. Moreover, commonly used model-based clustering algorithms do not generally provide a mechanism to quantify the relative contribution of each gene to the final partitioning of the data. We propose a penalised mixture of Student's t distributions for model-based clustering and gene ranking. Together with a bootstrap procedure, the proposed approach provides a means for ranking genes according to their contributions to the clustering process. Experimental results show that the algorithm performs well comparably to traditional Gaussian mixtures in the presence of outliers and longer tailed distributions. The algorithm also identifies the true informative genes with high sensitivity, and achieves improved model selection. An illustrative application to breast cancer data is also presented which confirms established tumor subclasses.

preprint2012arXiv

Static Parameter Estimation for ABC Approximations of Hidden Markov Models

In this article we focus on Maximum Likelihood estimation (MLE) for the static parameters of hidden Markov models (HMMs). We will consider the case where one cannot or does not want to compute the conditional likelihood density of the observation given the hidden state because of increased computational complexity or analytical intractability. Instead we will assume that one may obtain samples from this conditional likelihood and hence use approximate Bayesian computation (ABC) approximations of the original HMM. ABC approximations are biased, but the bias can be controlled to arbitrary precision via a parameter ε>0; the bias typically goes to zero as ε\searrow 0. We first establish that the bias in the log-likelihood and gradient of the log-likelihood of the ABC approximation, for a fixed batch of data, is no worse than \mathcal{O}(nε), n being the number of data; hence, for computational reasons, one might expect reasonable parameter estimates using such an ABC approximation. Turning to the computational problem of estimating $θ$, we propose, using the ABC-sequential Monte Carlo (SMC) algorithm in Jasra et al. (2012), an approach based upon simultaneous perturbation stochastic approximation (SPSA). Our method is investigated on two numerical examples

preprint2011arXiv

Error Bounds and Normalizing Constants for Sequential Monte Carlo in High Dimensions

In a recent paper Beskos et al (2011), the Sequential Monte Carlo (SMC) sampler introduced in Del Moral et al (2006), Neal (2001) has been shown to be asymptotically stable in the dimension of the state space d at a cost that is only polynomial in d, when N the number of Monte Carlo samples, is fixed. More precisely, it has been established that the effective sample size (ESS) of the ensuing (approximate) sample and the Monte Carlo error of fixed dimensional marginals will converge as $d$ grows, with a computational cost of $\mathcal{O}(Nd^2)$. In the present work, further results on SMC methods in high dimensions are provided as $d\to\infty$ and with $N$ fixed. We deduce an explicit bound on the Monte-Carlo error for estimates derived using the SMC sampler and the exact asymptotic relative $\mathbb{L}_2$-error of the estimate of the normalizing constant. We also establish marginal propagation of chaos properties of the algorithm. The accuracy in high-dimensions of some approximate SMC-based filtering schemes is also discussed.

preprint2011arXiv

On nonlinear Markov chain Monte Carlo

Let $\mathscr{P}(E)$ be the space of probability measures on a measurable space $(E,\mathcal{E})$. In this paper we introduce a class of nonlinear Markov chain Monte Carlo (MCMC) methods for simulating from a probability measure $π\in\mathscr{P}(E)$. Nonlinear Markov kernels (see [Feynman--Kac Formulae: Genealogical and Interacting Particle Systems with Applications (2004) Springer]) $K:\mathscr{P}(E)\times E\rightarrow\mathscr{P}(E)$ can be constructed to, in some sense, improve over MCMC methods. However, such nonlinear kernels cannot be simulated exactly, so approximations of the nonlinear kernels are constructed using auxiliary or potentially self-interacting chains. Several nonlinear kernels are presented and it is demonstrated that, under some conditions, the associated approximations exhibit a strong law of large numbers; our proof technique is via the Poisson equation and Foster--Lyapunov conditions. We investigate the performance of our approximations with some simulations.

preprint2011arXiv

Parameter Estimation for Hidden Markov Models with Intractable Likelihoods

Approximate Bayesian computation (ABC) is a popular technique for approximating likelihoods and is often used in parameter estimation when the likelihood functions are analytically intractable. Although the use of ABC is widespread in many fields, there has been little investigation of the theoretical properties of the resulting estimators. In this paper we give a theoretical analysis of the asymptotic properties of ABC based maximum likelihood parameter estimation for hidden Markov models. In particular, we derive results analogous to those of consistency and asymptotic normality for standard maximum likelihood estimation. We also discuss how Sequential Monte Carlo methods provide a natural method for implementing likelihood based ABC procedures.

preprint2010arXiv

Robust and Adaptive Algorithms for Online Portfolio Selection

We present an online approach to portfolio selection. The motivation is within the context of algorithmic trading, which demands fast and recursive updates of portfolio allocations, as new data arrives. In particular, we look at two online algorithms: Robust-Exponentially Weighted Least Squares (R-EWRLS) and a regularized Online minimum Variance algorithm (O-VAR). Our methods use simple ideas from signal processing and statistics, which are sometimes overlooked in the empirical financial literature. The two approaches are evaluated against benchmark allocation techniques using 4 real datasets. Our methods outperform the benchmark allocation techniques in these datasets, in terms of both computational demand and financial performance.

preprint2010arXiv

Sequential Monte Carlo Methods for Option Pricing

In the following paper we provide a review and development of sequential Monte Carlo (SMC) methods for option pricing. SMC are a class of Monte Carlo-based algorithms, that are designed to approximate expectations w.r.t a sequence of related probability measures. These approaches have been used, successfully, for a wide class of applications in engineering, statistics, physics and operations research. SMC methods are highly suited to many option pricing problems and sensitivity/Greek calculations due to the nature of the sequential simulation. However, it is seldom the case that such ideas are explicitly used in the option pricing literature. This article provides an up-to date review of SMC methods, which are appropriate for option pricing. In addition, it is illustrated how a number of existing approaches for option pricing can be enhanced via SMC. Specifically, when pricing the arithmetic Asian option w.r.t a complex stochastic volatility model, it is shown that SMC methods provide additional strategies to improve estimation.

preprint2010arXiv

The Time Machine: A Simulation Approach for Stochastic Trees

In the following paper we consider a simulation technique for stochastic trees. One of the most important areas in computational genetics is the calculation and subsequent maximization of the likelihood function associated to such models. This typically consists of using importance sampling (IS) and sequential Monte Carlo (SMC) techniques. The approach proceeds by simulating the tree, backward in time from observed data, to a most recent common ancestor (MRCA). However, in many cases, the computational time and variance of estimators are often too high to make standard approaches useful. In this paper we propose to stop the simulation, subsequently yielding biased estimates of the likelihood surface. The bias is investigated from a theoretical point of view. Results from simulation studies are also given to investigate the balance between loss of accuracy, saving in computing time and variance reduction.