When a hidden Markov model permits the conditional likelihood of an observation given the hidden process to be zero, all particle simulations from one observation time to the next could produce zeros. If so, the filtering distribution cannot be estimated and the estimated parameter likelihood is zero. The alive particle filter addresses this by simulating a random number of particles for each inter-observation interval, stopping after a target number of non-zero conditional likelihoods. For outlying observations or poor parameter values, a non-zero result can be extremely unlikely, and computational costs prohibitive. We introduce the Frankenfilter, a principled, partially alive particle filter that targets a user-defined amount of success whilst fixing lower and upper bounds on the number of simulations. The Frankenfilter produces unbiased estimators of the likelihood, suitable for pseudo-marginal Metropolis–Hastings (PMMH). We demonstrate that PMMH with the Frankenfilter is more robust to outliers and mis-specified initial parameter values than PMMH using standard particle filters, and is typically at least 2-3 times more efficient. We also provide advice for choosing the amount of success. In the case of n exact observations, this is particularly simple: target n successes.

The Metropolis-Hastings (MH) algorithm is one of the most widely used Markov Chain Monte Carlo schemes for generating samples from Bayesian posterior distributions. The algorithm is asymptotically exact, flexible and easy to implement. However, in the context of Bayesian inference for large datasets, evaluating the likelihood on the full data for thousands of iterations until convergence can be prohibitively expensive. This paper introduces a new subsample MH algorithm that satisfies detailed balance with respect to the target posterior and utilises control variates to enable exact, efficient Bayesian inference on datasets with large numbers of observations. Through theoretical results, simulation experiments and real-world applications on certain generalised linear models, we demonstrate that our method requires substantially smaller subsamples and is computationally more efficient than the standard MH algorithm and other exact subsample MH algorithms.

The ensemble Kalman filter (EnKF) is a popular technique for performing inference in state-space models (SSMs), particularly when the dynamic process is high-dimensional. Unlike reweighting methods such as sequential Monte Carlo (SMC, i.e. particle filters), the EnKF leverages either the linear Gaussian structure of the SSM or an approximation thereof, to maintain diversity of the sampled latent states (the so-called ensemble members) via shifting-based updates. Joint parameter and state inference using an EnKF is typically achieved by augmenting the state vector with the static parameter. In this case, it is assumed that both parameters and states follow a linear Gaussian state-space model, which may be unreasonable in practice. In this paper, we combine the reweighting and shifting methods by replacing the particle filter used in the SMC2 algorithm of Chopin et al. (2013), with the ensemble Kalman filter. Hence, parameter particles are weighted according to the estimated observed-data likelihood from the latest observation computed by the EnKF, and particle diversity is maintained via a resample-move step that targets the marginal parameter posterior under the EnKF. Extensions to the resulting algorithm are proposed, such as the use of a delayed acceptance kernel in the rejuvenation step and incorporation of nonlinear observation models. We illustrate the resulting methodology via several applications.

A state-space model is a time-series model that has an unobserved latent process from which we take noisy measurements over time. The observations are conditionally independent given the latent process and the latent process itself is Markovian. These properties lead to simplifications for the conditional distribution of the latent process given the parameters and the observations. This chapter looks at how we can leverage the properties of state-space models to construct efficient MCMC samplers. We consider a range of Gibbs-sampler schemes, including those which use the forward-backward algorithm to simulate from the full conditional of the latent process given the parameters. For models where the forward-backward algorithm is not applicable we look at particle MCMC algorithms that, given the parameters, use particle filters to approximately simulate from the latent process or estimate the likelihood of the observations. Throughout, we provide intuition and informally discuss theory about the properties of the model that impact the efficiency of the different algorithms and how approaches such as reparameterization can improve mixing.

A graduate-level introduction to advanced topics in Markov chain Monte Carlo (MCMC), as applied broadly in the Bayesian computational context. The topics covered have emerged as recently as the last decade and include stochastic gradient MCMC, non-reversible MCMC, continuous time MCMC, and new techniques for convergence assessment. A particular focus is on cutting-edge methods that are scalable with respect to either the amount of data, or the data dimension, motivated by the emerging high-priority application areas in machine learning and AI. Examples are woven throughout the text to demonstrate how scalable Bayesian learning methods can be implemented. This text could form the basis for a course and is sure to be an invaluable resource for researchers in the field.

In medical research, understanding changes in outcome measurements is crucial for inferring shifts in health conditions. However, traditional methods often struggle with large, irregularly longitudinal data and fail to account for the tendency of individuals in poorer health to interact more frequently with the healthcare system. Additionally, clinical data can lack information on terminating events like death. To address these challenges, we start from the continuous-time hidden Markov model which models observed data as outcomes influenced by latent health states. Our extension incorporates a point process to account for the impact of health states on observation timings and includes a "death" state to model unobserved terminating events through a Poisson process, where transition rates depend on the latent health state. This approach captures both the severity of the disease and the timing of healthcare interactions. We present an exact Gibbs sampler procedure that alternates between sampling the latent health state paths and the model parameters. By including the "death" state, we mitigate biases in parameter estimation that would arise from solely modelling "live" health states. Simulation studies demonstrate that the proposed Gibbs sampler performs effectively. We apply our method to Canadian healthcare data, offering valuable insights for healthcare management.

The k-year return levels of insurance losses due to flooding can be estimated by simulating and then summing a large number of independent losses for each of a large number of hypothetical years of flood events, and replicating this a large number of times. This leads to repeated realisations of the total losses over each year in a long sequence of years, from which return levels and their uncertainty can be estimated; the procedure, however, is highly computationally intensive. We develop and use a new, Bennett-like concentration inequality in a procedure that provides conservative but relatively accurate estimates of return levels at a fraction of the computational cost. Bennett's inequality provides concentration bounds on deviations of a sum of independent random variables from its expectation; it accounts for the different variances of each of the variables but uses only a single, uniform upper bound on their support. Motivated by the variability in the total insured value of insurance risks within a portfolio, we consider the case where the bounds on the support can vary by an order of magnitude or more, and obtain tractable concentration bounds. Simulation studies and application to a representative portfolio demonstrate the substantial improvement of our bounds over those obtained through Bennett's inequality. We then develop an importance sampling procedure that repeatedly samples the loss for each year from the distribution implied by the concentration inequality, leading to conservative estimates of the return levels and their uncertainty.

We congratulate the authors on an important contribution to the exciting area of generative models: providing a simple framework for applying the idea of diffusion generative models to data defined on general spaces and for general noising processes.Diffusion models have shown remarkable empirical performance for several applications, most notably as a way of defining implicit models for images.However it is unclear, to us at least, why these methods are so successful.They typically involve fitting highly-parameterised models for the score, and it is unclear why one would expect to be able to find good parameter estimates for such models.Furthermore, such models may be sufficiently flexible to learn the exact denoising process for the samples they are trained on-thereby collapsing the diffusion model to a categorical sampler from the training data.Could the authors provide insight into why these methods work?Is there intuition as to what types of data their denoising Markov models will work well for?Diffusion models are often used when the data are supported on a low-dimensional manifold of the state-space.However, the KL divergence which these models implicitly minimize can be infinite when the supports of the distributions are disjoint.Can the proposed framework extend to objectives which do not have this issue?Or is the use of KL beneficial, in that it actually helps force the model to sample from close to the manifold?The proposed framework parameterises the generative process in terms of a function β(x; t), whose optimal value is the marginal density of the noised data at time t.However, in practice the authors use different parameterisations, mainly based on differences in log β(x; t).Can the authors give more intuition as to why this is the best choice?And do they think this would be the case in general?Finally, both parameterisations of the continuous-time Markov chain model revolve around log r(x, y; t) = log β(x; t) -log β(y; t).As far as we can tell, the authors do not 1

The general applicability and ease of use of the pseudo-marginal Metropolis–Hastings (PMMH) algorithm, and particle Metropolis–Hastings in particular, makes it a popular method for inference on discretely observed Markovian stochastic processes. The performance of these algorithms and, in the case of particle Metropolis–Hastings, the trade off between improved mixing through increased accuracy of the estimator and the computational cost were investigated independently in two papers, both published in 2015. Each suggested choosing the number of particles so that the variance of the logarithm of the estimator of the posterior at a fixed sensible parameter value is approximately 1. This advice has been widely and successfully adopted. We provide new, remarkably simple upper and lower bounds on the asymptotic variance of PMMH algorithms. The bounds explain how blindly following the 2015 advice can hide serious issues with the algorithm and they strongly suggest an alternative criterion. In most situations our guidelines and those from 2015 closely coincide; however, when the two differ it is safer to follow the new guidance. An extension of one of our bounds shows how the use of correlated proposals can fundamentally shift the properties of pseudo-marginal algorithms, so that asymptotic variances that were infinite under the PMMH kernel become finite.

There has been a recent surge of interest in coupling methods for Markov chain Monte Carlo algorithms: they facilitate convergence quantification and unbiased estimation, while exploiting embarrassingly parallel computing capabilities. Motivated by these, we consider the design and analysis of couplings of the random walk Metropolis algorithm which scale well with the dimension of the target measure. Methodologically, we introduce a low-rank modification of the synchronous coupling that is provably optimally contractive in standard high-dimensional asymptotic regimes. We expose a shortcoming of the reflection coupling, the state of the art at the time of writing, and we propose a modification which mitigates the issue. Our analysis bridges the gap to the optimal scaling literature and builds a framework of asymptotic optimality which may be of independent interest. We illustrate the applicability of our proposed couplings, and the potential for extending our ideas, with various numerical experiments.

This book aims to provide a graduate-level introduction to advanced topics in Markov chain Monte Carlo (MCMC) algorithms, as applied broadly in the Bayesian computational context. Most, if not all of these topics (stochastic gradient MCMC, non-reversible MCMC, continuous time MCMC, and new techniques for convergence assessment) have emerged as recently as the last decade, and have driven substantial recent practical and theoretical advances in the field. A particular focus is on methods that are scalable with respect to either the amount of data, or the data dimension, motivated by the emerging high-priority application areas in machine learning and AI.

We present new methodologies for Bayesian inference on the rate parameters of a discretely observed continuous-time Markov jump processes with a countably infinite state space. The usual method of choice for inference, particle Markov chain Monte Carlo (particle MCMC), struggles when the observation noise is small. We consider the most challenging regime of exact observations and provide two new methodologies for inference in this case: the minimal extended state space algorithm (MESA) and the nearly minimal extended state space algorithm (nMESA). By extending the Markov chain Monte Carlo state space, both MESA and nMESA use the exponentiation of finite rate matrices to perform exact Bayesian inference on the Markov jump process even though its state space is countably infinite. Numerical experiments show improvements over particle MCMC of between a factor of three and several orders of magnitude.

Divide-and-conquer strategies for Monte Carlo algorithms are an increasingly popular approach to making Bayesian inference scalable to large data sets. In its simplest form, the data are partitioned across multiple computing cores and a separate Markov chain Monte Carlo algorithm on each core targets the associated partial posterior distribution, which we refer to as a sub-posterior, that is, the posterior given only the data from the segment of the partition associated with that core. Divide-and-conquer techniques reduce computational, memory and disk bottle necks but make it difficult to recombine the sub-posterior samples. We propose SwISS: Sub-posteriors with Inflation, Scaling and Shifting; a new approach for recombining the sub-posterior samples, which is simple to apply, scales to high-dimensional parameter spaces and accurately approximates the original posterior distribution through affine transformations of the sub-posterior samples. We prove that our transformation is asymptotically optimal across a natural set of affine transformations and illustrate the efficacy of SwISS against competing algorithms on synthetic and real-world data sets.

We introduced the Hug and Hop Markov chain Monte Carlo algorithm for estimating expectations with respect to an intractable distribution. The algorithm alternates between two kernels: Hug and Hop. Hug is a non-reversible kernel that repeatedly applies the bounce mechanism from the recently proposed Bouncy Particle Sampler to produce a proposal point far from the current position, yet on almost the same contour of the target density, leading to a high acceptance probability. Hug is complemented by Hop, which deliberately proposes jumps between contours and has an efficiency that degrades very slowly with increasing dimension. There are many parallels between Hug and Hamiltonian Monte Carlo using a leapfrog integrator, including the order of the integration scheme, however Hug is also able to make use of local Hessian information without requiring implicit numerical integration steps, and its performance is not terminally affected by unbounded gradients of the log-posterior. We test Hug and Hop empirically on a variety of toy targets and real statistical models and find that it can, and often does, outperform Hamiltonian Monte Carlo.

Stochastic kinetic models (SKMs) are increasingly used to account for the inherent stochasticity exhibited by interacting populations of species in areas such as epidemiology, population ecology and systems biology. Species numbers are modelled using a continuous-time stochastic process, and, depending on the application area of interest, this will typically take the form of a Markov jump process or an Itô diffusion process. Widespread use of these models is typically precluded by their computational complexity. In particular, performing exact fully Bayesian inference in either modelling framework is challenging due to the intractability of the observed data likelihood, necessitating the use of computationally intensive techniques such as particle Markov chain Monte Carlo (particle MCMC). It is proposed to increase the computational and statistical efficiency of this approach by leveraging the tractability of an inexpensive surrogate derived directly from either the jump or diffusion process. The surrogate is used in three ways: in the design of a gradient-based parameter proposal, to construct an appropriate bridge and in the first stage of a delayed-acceptance step. The resulting approach, which exactly targets the posterior of interest, offers substantial gains in efficiency over a standard particle MCMC implementation.

Amongst Markov chain Monte Carlo algorithms, Hamiltonian Monte Carlo (HMC) is often the algorithm of choice for complex, high-dimensional target distributions; however, its efficiency is notoriously sensitive to the choice of the integration-time tuning parameter, $T$. When integrating both forward and backward in time using the same leapfrog integration step as HMC, the set of local maxima in the potential along a path, or apogees, is the same whatever point (position and momentum) along the path is chosen to initialise the integration. We present the Apogee to Apogee Path Sampler (AAPS), which utilises this invariance to create a simple yet generic methodology for constructing a path, proposing a point from it and accepting or rejecting that proposal so as to target the intended distribution. We demonstrate empirically that AAPS has a similar efficiency to HMC but is much more robust to the setting of its equivalent tuning parameter, the number of apogees that the path crosses.

Most Markov chain Monte Carlo methods operate in discrete time and are reversible with respect to the target probability. Nevertheless, it is now understood that the use of nonreversible Markov chains can be beneficial in many contexts. In particular, the recently proposed bouncy particle sampler leverages a continuous-time and nonreversible Markov process, and empirically shows state-of-the-art performance when used to explore certain probability densities; however, its implementation typically requires the computation of local upper bounds on the gradient of the log target density. We present the discrete bouncy particle sampler, a general algorithm based on a guided random walk, a partial refreshment of direction and a delayed-rejection step. We show that the bouncy particle sampler can be understood as a scaling limit of a special case of our algorithm. In contrast to the bouncy particle sampler, implementing the discrete bouncy particle sampler only requires pointwise evaluation of the target density and its gradient. We propose extensions of the basic algorithm for situations when the exact gradient of the target density is not available. In a Gaussian setting, we establish a scaling limit for the radial process as the dimension increases to infinity. We leverage this result to obtain the theoretical efficiency of the discrete bouncy particle sampler as a function of the partial-refreshment parameter, which leads to a simple and robust tuning criterion. A further analysis in a more general setting suggests that this tuning criterion applies more generally. Theoretical and empirical efficiency curves are then compared for different targets and algorithm variations.

Summary We introduce a general framework for monitoring, modeling, and predicting the recruitment to multi-center clinical trials. The work is motivated by overly optimistic and narrow prediction intervals produced by existing time-homogeneous recruitment models for multi-center recruitment. We first present two tests for detection of decay in recruitment rates, together with a power study. We then introduce a model based on the inhomogeneous Poisson process with monotonically decaying intensity, motivated by recruitment trends observed in oncology trials. The general form of the model permits adaptation to any parametric curve-shape. A general method for constructing sensible parameter priors is provided and Bayesian model averaging is used for making predictions which account for the uncertainty in both the parameters and the model. The validity of the method and its robustness to misspecification are tested using simulated datasets. The new methodology is then applied to oncology trial data, where we make interim accrual predictions, comparing them to those obtained by existing methods, and indicate where unexpected changes in the accrual pattern occur.

The plug-in estimator of the Wasserstein distance is known to be conservative, however its usefulness is severely limited when the distributions are similar as its bias does not decay to zero with the true Wasserstein distance. We propose a linear combination of plug-in estimators for the squared 2-Wasserstein distance with a reduced bias that decays to zero with the true distance. The new estimator is provably conservative provided one distribution is appropriately overdispersed with respect the other, and is unbiased when the distributions are equal. We apply it to approximately bound from above the 2-Wasserstein distance between the target and current distribution in Markov chain Monte Carlo, running multiple identically distributed chains which start, and remain, overdispersed with respect to the target. Our bound consistently outperforms the current state-of-the-art bound, which uses coupling, improving mixing time bounds by up to an order of magnitude.

A delayed-acceptance version of a Metropolis–Hastings algorithm can be useful for Bayesian inference when it is computationally expensive to calculate the true posterior, but a computationally cheap approximation is available; the delayed-acceptance kernel targets the same posterior as its associated “parent” Metropolis-Hastings kernel. Although the asymptotic variance of the ergodic average of any functional of the delayed-acceptance chain cannot be less than that obtained using its parent, the average computational time per iteration can be much smaller and so for a given computational budget the delayed-acceptance kernel can be more efficient. When the asymptotic variance of the ergodic averages of all $$L^2$$ functionals of the chain are finite, the kernel is said to be variance bounding. It has recently been noted that a delayed-acceptance kernel need not be variance bounding even when its parent is. We provide sufficient conditions for inheritance: for non-local algorithms, such as the independence sampler, the discrepancy between the log density of the approximation and that of the truth should be bounded; for local algorithms, two alternative sets of conditions are provided. As a by-product of our initial, general result we also supply sufficient conditions on any pair of proposals such that, for any shared target distribution, if a Metropolis-Hastings kernel using one of the proposals is variance bounding then so is the Metropolis-Hastings kernel using the other proposal.