Consider a population evolving as a discrete-time supercritical multi-type Galton–Watson process. Suppose we run the process for T generations, then sample k individuals uniformly at generation T and trace their genealogy backwards in time. In the limiting regime as T →∞, the expected behaviour of the sample's ancestry has been analysed extensively in the single-type case and, more recently, for multi-type processes in the critical case. In this paper, we present a formula for the distribution function of the generation t of the most recent common ancestor in terms of the limiting distribution of the normalised population size. In addition, we provide effective bounds for the decay of this distribution function to 1 in terms of the harmonic moments of the population size at generation t. In order to better understand the behaviour of these harmonic moments, we use a multi-type generalisation of the Harris–Sevastyanov transformation to express harmonic moments at generation t in terms of moments of the transformed process at the first generation. We present numerical results demonstrating that it is possible to approximate the coalescence time distribution effectively in practical settings.

Many real-world processes can be represented as compositions of functions along a directed acyclic graph (DAG). In causal modelling, these correspond to the underlying mechanisms; in engineering, to multiple fidelity levels; and in gene-regulatory networks, to transcription factors. These functions are partially observed across the DAG, with noisy and heterogeneously sampled measurements, posing significant challenges for reconstruction, uncertainty propagation, and inference. To tackle these challenges, we place priors over functions and naturally arrive at Deep Gaussian Processes over DAGs. We theoretically study their prior-collapse behaviour, and the effect of graph topology and intermediate observations on the preservation of information. We obtain almost-sure lower bounds on the asymptotic frequency of depths at which the distinction between inputs is preserved, identify broad kernel classes for which these hold, and prove an observation by on the role of input connections. We offer a structured variational approximation that retains graph dependencies, preserves compositional uncertainty, and captures the explaining-away behaviour of colliders. Finally, we empirically validate our theoretical results and our methodology, and model a latent-collider DAG, a protein signalling network, and a multi-fidelity heavy-ion collision emulation task, attaining state-of-the-art performance while recovering low-fidelity contributions and yielding interpretability of the simulator hierarchy.

Conditional sequential Monte Carlo (CSMC) algorithms arise in particle Markov chain Monte Carlo and a number of related settings. As in standard sequential Monte Carlo (SMC) algorithms, it is possible to employ a number of approaches to resampling within CSMC, but some additional care is required to arrive at a valid algorithm. We present a simple framework for implementing valid SMC and CSMC algorithms which (a) covers most known resampling schemes including those with a complicated dependence structure like systematic resampling, but also adaptive resampling, and even more `exotic' schemes like a version of chopthin resampling; (b) explains how to implement conditional analogues of these and other well known resampling schemes without randomly permuting/shifting the order of the ancestor indices; (c) requires only very weak assumptions which include neither (marginal) "unbiasedness" nor exchangeability.

We show that genealogical trees arising from a broad class of non-neutral models of population evolution converge to the Kingman coalescent under a suitable rescaling of time. As well as non-neutral biological evolution, our results apply to genetic algorithms encompassing the prominent class of sequential Monte Carlo (SMC) methods. The time rescaling we need differs slightly from that used in classical results for convergence to the Kingman coalescent, which has implications for the performance of different resampling schemes in SMC algorithms. In addition, our work substantially simplifies earlier proofs of convergence to the Kingman coalescent, and corrects an error common to several earlier results.

We introduce methodology for real-time inference in general-state-space hidden Markov models. Specifically, we extend recent advances in controlled sequential Monte Carlo (CSMC) methods-originally proposed for offline smoothing-to the online setting via a rolling window mechanism. Our novel online rolling controlled sequential Monte Carlo (ORCSMC) algorithm employs two particle systems to simultaneously estimate twisting functions and perform filtering, ensuring real-time adaptivity to new observations while maintaining bounded computational cost. Numerical results on linear-Gaussian, stochastic volatility, and neuroscience models demonstrate improved estimation accuracy and robustness in higher dimensions, compared to standard particle filtering approaches. The method offers a statistically efficient and practical solution for sequential and real-time inference in complex latent variable models.

We present a new framework for analysing the expectation-maximization (em) algorithm. Drawing on recent advances in the theory of gradient flows over Euclidean-Wasserstein spaces, we extend techniques from alternating minimization in Euclidean spaces to the em algorithm, via its representation as coordinatewise minimization of the free energy. In so doing, we obtain finite-sample error bounds and exponential convergence of the em algorithm under a natural generalization of the log-Sobolev inequality. We further show that this framework naturally extends to several variants of the em algorithm, offering a unified approach for studying such algorithms.

We prove non-asymptotic error bounds for particle gradient descent (PGD)~(Kuntz et al., 2023), a recently introduced algorithm for maximum likelihood estimation of large latent variable models obtained by discretizing a gradient flow of the free energy. We begin by showing that, for models satisfying a condition generalizing both the log-Sobolev and the Polyak--{\L}ojasiewicz inequalities (LSI and P{\L}I, respectively), the flow converges exponentially fast to the set of minimizers of the free energy. We achieve this by extending a result well-known in the optimal transport literature (that the LSI implies the Talagrand inequality) and its counterpart in the optimization literature (that the P{\L}I implies the so-called quadratic growth condition), and applying it to our new setting. We also generalize the Bakry--\'Emery Theorem and show that the LSI/P{\L}I generalization holds for models with strongly concave log-likelihoods. For such models, we further control PGD's discretization error, obtaining non-asymptotic error bounds. While we are motivated by the study of PGD, we believe that the inequalities and results we extend may be of independent interest.

We introduce an inferential framework for a wide class of semi-linear stochastic differential equations (SDEs). Recent work has shown that numerical splitting schemes can preserve critical properties of such types of SDEs, give rise to explicit pseudolikelihoods, and hence allow for parameter inference for fully observed processes. Here, under several discrete time observation regimes (particularly, partially and fully observed with and without noise), we represent the implied pseudolikelihood as the normalising constant of a Feynman–Kac flow, allowing its efficient estimation via controlled sequential Monte Carlo and adapt likelihood-based methods to exploit this pseudolikelihood for inference. The strategy developed herein allows us to obtain good inferential results across a range of problems. Using diffusion bridges, we are able to computationally reduce bias coming from time-discretisation without recourse to more complex numerical schemes which typically require considerable application-specific efforts. Simulations illustrate that our method provides an excellent trade-off between computational efficiency and accuracy, under hypoellipticity, for both point and posterior estimation. Application to a neuroscience example shows the good performance of the method in challenging settings.

We introduce a new class of sequential Monte Carlo methods which reformulates the essence of the nested sampling (NS) method of Skilling in terms of sequential Monte Carlo techniques. Two new algorithms are proposed: nested sampling via sequential Monte Carlo (NS-SMC) and adaptive nested sampling via sequential Monte Carlo (ANS-SMC). The new framework allows convergence results to be obtained in the setting when Markov chain Monte Carlo (MCMC) is used to produce new samples. An additional benefit is that marginal-likelihood (normalizing constant) estimates given by NS-SMC are unbiased. In contrast to NS, the analysis of our proposed algorithms does not require the (unrealistic) assumption that the simulated samples be independent. We show that a minor adjustment to our ANS-SMC algorithm recovers the original NS algorithm, which provides insights as to why NS seems to produce accurate estimates despite a typical violation of its assumptions. A numerical study is conducted where the performance of the proposed algorithms and temperature-annealed SMC is compared on challenging problems. Code for the experiments is made available online at https://github.com/LeahPrice/SMC-NS.

Maximum likelihood estimation (MLE) of latent variable models is often recast as an optimization problem over the extended space of parameters and probability distributions. For example, the Expectation Maximization (EM) algorithm can be interpreted as coordinate descent applied to a suitable free energy functional over this space. Recently, this perspective has been combined with insights from optimal transport and Wasserstein gradient flows to develop particle-based algorithms applicable to wider classes of models than standard EM. Drawing inspiration from prior works which interpret `momentum-enriched' optimisation algorithms as discretizations of ordinary differential equations, we propose an analogous dynamical systems-inspired approach to minimizing the free energy functional over the extended space of parameters and probability distributions. The result is a dynamic system that blends elements of Nesterov's Accelerated Gradient method, the underdamped Langevin diffusion, and particle methods. Under suitable assumptions, we establish quantitative convergence of the proposed system to the unique minimiser of the functional in continuous time. We then propose a numerical discretization of this system which enables its application to parameter estimation in latent variable models. Through numerical experiments, we demonstrate that the resulting algorithm converges faster than existing methods and compares favourably with other (approximate) MLE algorithms.

Semi-implicit variational inference (SIVI) enriches the expressiveness of variational families by utilizing a kernel and a mixing distribution to hierarchically define the variational distribution. Existing SIVI methods parameterize the mixing distribution using implicit distributions, leading to intractable variational densities. As a result, directly maximizing the evidence lower bound (ELBO) is not possible, so they resort to one of the following: optimizing bounds on the ELBO, employing costly inner-loop Markov chain Monte Carlo runs, or solving minimax objectives. In this paper, we propose a novel method for SIVI called Particle Variational Inference (PVI) which employs empirical measures to approximate the optimal mixing distributions characterized as the minimizer of a free energy functional. PVI arises naturally as a particle approximation of a Euclidean–Wasserstein gradient flow and, unlike prior works, it directly optimizes the ELBO whilst making no parametric assumption about the mixing distribution. Our empirical results demonstrate that PVI performs favourably compared to other SIVI methods across various tasks. Moreover, we provide a theoretical analysis of the behaviour of the gradient flow of a related free energy functional: establishing the existence and uniqueness of solutions as well as propagation of chaos results.

Solving Fredholm equations of the first kind is crucial in many areas of the applied sciences. In this work we consider integral equations featuring kernels which may be expressed as scalar multiples of conservative (i.e. Markov) kernels and we adopt a variational point of view by considering a minimization problem in the space of probability measures with an entropic regularization. Contrary to classical approaches which discretize the domain of the solutions, we introduce an algorithm to asymptotically sample from the unique solution of the regularized minimization problem. As a result our estimators do not depend on any underlying grid and have better scalability properties than most existing methods. Our algorithm is based on a particle approximation of the solution of a McKean–Vlasov stochastic differential equation associated with the Wasserstein gradient flow of our variational formulation. We prove the convergence towards a minimizer and provide practical guidelines for its numerical implementation. Finally, our method is compared with other approaches on several examples including density deconvolution and epidemiology.

The problem of incorporating information from observations received serially in time is widespread in the field of uncertainty quantification. Within a probabilistic framework, such problems can be addressed using standard filtering techniques. However, in many real-world problems, some (or all) of the uncertainty is epistemic, arising from a lack of knowledge, and is difficult to model probabilistically. This paper introduces a possibilistic ensemble Kalman filter designed for this setting and characterizes some of its properties. Using possibility theory to describe epistemic uncertainty is appealing from a philosophical perspective, and it is easy to justify certain heuristics often employed in standard ensemble Kalman filters as principled approaches to capturing uncertainty within it. The possibilistic approach motivates a robust mechanism for characterizing uncertainty which shows good performance with small sample sizes, and can outperform standard ensemble Kalman filters at given sample size, even when dealing with genuinely aleatoric uncertainty.

We propose a divide-and-conquer approach to filtering which decomposes the state variable into low-dimensional components to which standard particle filtering tools can be successfully applied and recursively merges them to recover the full filtering distribution. It is less dependent upon factorization of transition densities and observation likelihoods than competing approaches and can be applied to a broader class of models. Performance is compared with state-of-the-art methods on a benchmark problem and it is demonstrated that the proposed method is broadly comparable in settings in which those methods are applicable, and that it can be applied in settings in which they cannot.

We provide a comprehensive characterisation of the theoretical properties of the divide-and-conquer sequential Monte Carlo (DaC-SMC) algorithm. We firmly establish it as a well-founded method by showing that it possesses the same basic properties as conventional sequential Monte Carlo (SMC) algorithms do. In particular, we derive pertinent laws of large numbers, $L^p$ inequalities, and central limit theorems; and we characterize the bias in the normalized estimates produced by the algorithm and argue the absence thereof in the unnormalized ones. We further consider its practical implementation and several interesting variants; obtain expressions for its globally and locally optimal intermediate targets, auxiliary measures, and proposal kernels; and show that, in comparable conditions, DaC-SMC proves more statistically efficient than its direct SMC analogue. We close the paper with a discussion of our results, open questions, and future research directions.

In Bayesian inverse problems, one aims at characterizing the posterior distribution of a set of unknowns, given indirect measurements. For non-linear/non-Gaussian problems, analytic solutions are seldom available: Sequential Monte Carlo samplers offer a powerful tool for approximating complex posteriors, by constructing an auxiliary sequence of densities that smoothly reaches the posterior. Often the posterior depends on a scalar hyper-parameter, for which limited prior information is available. In this work, we show that properly designed Sequential Monte Carlo (SMC) samplers naturally provide an approximation of the marginal likelihood associated with this hyper-parameter for free, i.e. at a negligible additional computational cost. The proposed method proceeds by constructing the auxiliary sequence of distributions in such a way that each of them can be interpreted as a posterior distribution corresponding to a different value of the hyper-parameter. This can be exploited to perform selection of the hyper-parameter in Empirical Bayes (EB) approaches, as well as averaging across values of the hyper-parameter according to some hyper-prior distribution in Fully Bayesian (FB) approaches. For FB approaches, the proposed method has the further benefit of allowing prior sensitivity analysis at a negligible computational cost. In addition, the proposed method exploits particles at all the (relevant) iterations, thus alleviating one of the known limitations of SMC samplers, i.e. the fact that all samples at intermediate iterations are typically discarded. We show numerical results for two distinct cases where the hyper-parameter affects only the likelihood: a toy example, where an SMC sampler is used to approximate the full posterior distribution; and a brain imaging example, where a Rao-Blackwellized SMC sampler is used to approximate the posterior distribution of a subset of parameters in a conditionally linear Gaussian model.

Fredholm integral equations of the first kind are the prototypical example of ill-posed linear inverse problems. They model, among other things, reconstruction of distorted noisy observations and indirect density estimation and also appear in instrumental variable regression. However, their numerical solution remains a challenging problem. Many techniques currently available require a preliminary discretization of the domain of the solution and make strong assumptions about its regularity. For example, the popular expectation maximization smoothing (EMS) scheme requires the assumption of piecewise constant solutions which is inappropriate for most applications. We propose here a novel particle method that circumvents these two issues. This algorithm can be thought of as a Monte Carlo approximation of the EMS scheme which not only performs an adaptive stochastic discretization of the domain but also results in smooth approximate solutions. We analyze the theoretical properties of the EMS iteration and of the corresponding particle algorithm. Compared to standard EMS, we show experimentally that our novel particle method provides state-of-the-art performance for realistic systems, including motion deblurring and reconstruction of cross-section images of the brain from positron emission tomography.

We provide a framework which admits a number of ``marginal'' sequential Monte Carlo (SMC) algorithms as particular cases -- including the marginal particle filter [Klaas et al., 2005, in: Proceedings of Uncertainty in Artificial Intelligence, pp. 308--315], , the independent particle filter [Lin et al., 2005, Journal of the American Statistical Association 100, pp. 1412--1421] and linear-cost Approximate Bayesian Computation SMC [Sisson et al., 2007, Proceedings of the National Academy of Sciences (USA) 104, pp. 1760--1765.]. We provide conditions under which such algorithms obey laws of large numbers and central limit theorems and provide some further asymptotic characterizations. Finally, it is shown that the asymptotic variance of a class of estimators associated with certain marginal SMC algorithms is never greater than that of the estimators provided by a standard SMC algorithm using the same proposal distributions.

(Neal and Hinton, 1998) recast maximum likelihood estimation of any given latent variable model as the minimization of a free energy functional $F$, and the EM algorithm as coordinate descent applied to $F$. Here, we explore alternative ways to optimize the functional. In particular, we identify various gradient flows associated with $F$ and show that their limits coincide with $F$'s stationary points. By discretizing the flows, we obtain practical particle-based algorithms for maximum likelihood estimation in broad classes of latent variable models. The novel algorithms scale to high-dimensional settings and perform well in numerical experiments.

Markov chain Monte Carlo (MCMC) algorithms are based on the construction of a Markov Chain with transition probabilities P_μ(x,·), where μ indicates an invariant distribution of interest. In this work, we look at these transition probabilities as functions of their invariant distributions, and we develop a notion of derivative in the invariant distribution of a MCMC kernel. We build around this concept a set of tools that we refer to as Markov Chain Monte Carlo Calculus. This allows us to compare Markov chains with different invariant distributions within a suitable class via what we refer to as mean value inequalities. We explain how MCMC Calculus provides a natural framework to study algorithms using an approximation of an invariant distribution, also illustrating how it suggests practical guidelines for MCMC algorithms efficiency. We conclude this work by showing how the tools developed can be applied to prove convergence of interacting and sequential MCMC algorithms, which arise in the context of particle filtering.