Tempering is a popular tool in Bayesian computation, being used to transform a posterior distribution p_1 into a reference distribution p_0 that is more easily approximated. Several algorithms exist that start by approximating p_0 and proceed through a sequence of intermediate distributions p_t until an approximation to p_1 is obtained. Our contribution reveals that high-quality approximation of terms up to p_1 is not essential, as knowledge of the intermediate distributions enables posterior quantities of interest to be extrapolated. Specifically, we establish conditions under which posterior expectations are determined by their associated tempered expectations on any non-empty t interval. Harnessing this result, we propose novel methodology for approximating posterior expectations based on extrapolation and smoothing of tempered expectations, which we implement as a post-processing variance-reduction tool for sequential Monte Carlo.

Diffusion models are typically trained using score matching, a learning objective agnostic to the underlying noising process that guides the model. This paper argues that Markov noising processes enjoy an advantage over alternatives, as the Markov operators that govern the noising process are well-understood. Specifically, by leveraging the spectral decomposition of the infinitesimal generator of the Markov noising process, we obtain parametric estimates of the score functions simultaneously for all marginal distributions, using only sample averages with respect to the data distribution. The resulting operator-informed score matching provides both a standalone approach to sample generation for low-dimensional distributions, as well as a recipe for better informed neural score estimators in high-dimensional settings.

The use of heuristics to assess the convergence and compress the output of Markov chain Monte Carlo can be sub-optimal in terms of the empirical approximations that are produced. Typically a number of the initial states are attributed to 'burn in' and removed, while the remainder of the chain is 'thinned' if compression is also required. In this paper, we consider the problem of retrospectively selecting a subset of states, of fixed cardinality, from the sample path such that the approximation provided by their empirical distribution is close to optimal. A novel method is proposed, based on greedy minimisation of a kernel Stein discrepancy, that is suitable when the gradient of the log-target can be evaluated and approximation using a small number of states is required. Theoretical results guarantee consistency of the method and its effectiveness is demonstrated in the challenging context of parameter inference for ordinary differential equations. Software is available in the Stein Thinning package in Python, R and MATLAB.

The use of heuristics to assess the convergence and compress the output of Markov chain Monte Carlo can be sub-optimal in terms of the empirical approximations that are produced. Typically a number of the initial states are attributed to ‘burn in’ and removed, while the remainder of the chain is ‘thinned’ if compression is also required. In this paper, we consider the problem of retrospectively selecting a subset of states, of fixed cardinality, from the sample path such that the approximation provided by their empirical distribution is close to optimal. A novel method is proposed, based on greedy minimisation of a kernel Stein discrepancy, that is suitable when the gradient of the log-target can be evaluated and approximation using a small number of states is required. Theoretical results guarantee consistency of the method and its effectiveness is demonstrated in the challenging context of parameter inference for ordinary differential equations. Software is available in the Stein Thinning package in Python , R and MATLAB .

Markov chain Monte Carlo is the engine of modern Bayesian statistics, being used to approximate the posterior and derived quantities of interest. Despite this, the issue of how the output from a Markov chain is post-processed and reported is often overlooked. Convergence diagnostics can be used to control bias via burn-in removal, but these do not account for (common) situations where a limited computational budget engenders a bias-variance trade-off. The aim of this article is to review state-of-the-art techniques for post-processing Markov chain output. Our review covers methods based on discrepancy minimisation, which directly address the bias-variance trade-off, as well as general-purpose control variate methods for approximating expected quantities of interest.

Several researchers have proposed minimisation of maximum mean discrepancy (MMD) as a method to quantise probability measures, i.e., to approximate a distribution by a representative point set. We consider sequential algorithms that greedily minimise MMD over a discrete candidate set. We propose a novel non-myopic algorithm and, in order to both improve statistical efficiency and reduce computational cost, we investigate a variant that applies this technique to a mini-batch of the candidate set at each iteration. When the candidate points are sampled from the target, the consistency of these new algorithms-and their mini-batch variants-is established. We demonstrate the algorithms on a range of important computational problems, including optimisation of nodes in Bayesian cubature and the thinning of Markov chain output.

Uncertainty quantification (UQ) is a vital step in using mathematical models and simulations to take decisions. The field of cardiac simulation has begun to explore and adopt UQ methods to characterize uncertainty in model inputs and how that propagates through to outputs or predictions; examples of this can be seen in the papers of this issue. In this review and perspective piece, we draw attention to an important and under-addressed source of uncertainty in our predictions-that of uncertainty in the model structure or the equations themselves. The difference between imperfect models and reality is termed model discrepancy, and we are often uncertain as to the size and consequences of this discrepancy. Here, we provide two examples of the consequences of discrepancy when calibrating models at the ion channel and action potential scales. Furthermore, we attempt to account for this discrepancy when calibrating and validating an ion channel model using different methods, based on modelling the discrepancy using Gaussian processes and autoregressive-moving-average models, then highlight the advantages and shortcomings of each approach. Finally, suggestions and lines of enquiry for future work are provided. This article is part of the theme issue 'Uncertainty quantification in cardiac and cardiovascular modelling and simulation'.

Patient-specific cardiac models are now being used to guide therapies. The increased use of patient-specific cardiac simulations in clinical care will give rise to the development of virtual cohorts of cardiac models. These cohorts will allow cardiac simulations to capture and quantify inter-patient variability. However, the development of virtual cohorts of cardiac models will require the transformation of cardiac modelling from small numbers of bespoke models to robust and rapid workflows that can create large numbers of models. In this review, we describe the state of the art in virtual cohorts of cardiac models, the process of creating virtual cohorts of cardiac models, and how to generate the individual cohort member models, followed by a discussion of the potential and future applications of virtual cohorts of cardiac models.This article is part of the theme issue ‘Uncertainty quantification in cardiac and cardiovascular modelling and simulation’.

The results of a series of theoretical studies are reported, examining the convergence rate for different approximate representations of α-stable distributions. Although they play a key role in modelling random processes with jumps and discontinuities, the use of α-stable distributions in inference often leads to analytically intractable problems. The LePage series, which is a probabilistic representation employed in this work, is used to transform an intractable, infinite-dimensional inference problem into a finite-dimensional (conditionally Gaussian) parametric problem. A major component of our approach is the approximation of the tail of this series by a Gaussian random variable. Standard statistical techniques, such as ExpectationMaximization (EM), Markov chain Monte Carlo, and Particle Filtering, can then be readily applied. In addition to the asymptotic normality of the tail of this series, we establish explicit, nonasymptotic bounds on the approximation error. Their proofs follow classical Fourier-analytic arguments, using Esséen's smoothing lemma. Specifically, we consider the distance between the distributions of: (i) the tail of the series and an appropriate Gaussian; (ii) the full series and the truncated series; and (iii) the full series and the truncated series with an added Gaussian term. In all three cases, sharp bounds are established, and the theoretical results are compared with the actual distances (computed numerically) in specific examples of symmetric αstable distributions. This analysis facilitates the selection of appropriate truncations in practice and offers theoretical guarantees for the accuracy of resulting estimates. One of the main conclusions obtained is that, for the purposes of inference, the use of a truncated series together with an approximately Gaussian error term has superior statistical properties and is likely a preferable choice in practice.

In this paper we introduce a new class of state space models based on shot-noise simulation representations of nonGaussian Lévy-driven linear systems, represented as stochastic differential equations. In particular a conditionally Gaussian version of the models is proposed that is able to capture heavy-tailed non-Gaussianity while retaining tractability for inference procedures. We focus on a canonical class of such processes, the α-stable Lévy processes, which retain important properties such as self-similarity and heavy-tails, while emphasizing that broader classes of non-Gaussian Lévy processes may be handled by similar methodology. An important feature is that we are able to marginalise both the skewness and the scale parameters of these challenging models from posterior probability distributions. The models are posed in continuous time and so are able to deal with irregular data arrival times. Example modelling and inference procedures are provided using Rao-Blackwellised sequential Monte Carlo applied to a two-dimensional Langevin model, and this is tested on real exchange rate data.

We report the results of several theoretical studies into the convergence rate for certain random series representations of α-stable random variables, which are motivated by and find application in modelling heavy-tailed noise in time series analysis, inference, and stochastic processes. The use of α-stable noise distributions generally leads to analytically intractable inference problems. The particular version of the Poisson series representation invoked here implies that the resulting distributions are “conditionally Gaussian,” for which inference is relatively straightforward, although an infinite series is still involved. Our approach is to approximate the residual (or “tail”) part of the series from some point, c > 0, say, to∞, as a Gaussian random variable. Empirically, this approximation has been found to be very accurate for large c. We study the rate of convergence, as c → ∞, of this Gaussian approximation. This allows the selection of appropriate truncation parameters, so that a desired level of accuracy for the approximate model can be achieved. Explicit, nonasymptotic bounds are obtained for the Kolmogorov distance between the relevant distribution functions, through the application of probability-theoretic tools. The theoretical results obtained are found to be in very close agreement with numerical results obtained in earlier work.

We report the results of a series of numerical studies examining the convergence rate for some approximate representations of α-stable distributions, which are a highly intractable class of distributions for inference purposes. Our proposed representation turns the intractable inference for an infinite-dimensional series of parameters into an (approximately) conditionally Gaussian representation, to which standard inference procedures such as Expectation-Maximization (EM), Markov chain Monte Carlo (MCMC) and Particle Filtering can be readily applied. While we have previously proved the asymptotic convergence of this representation, here we study the rate of this convergence for finite values of a truncation parameter, c. This allows the selection of appropriate truncations for different parameter configurations and for the accuracy required for the model. The convergence is examined directly in terms of cumulative distribution functions and densities, through the application of the Berry theorems and Parseval theorems. Our results indicate that the behaviour of our representations is significantly superior to that of representations that simply truncate the series with no Gaussian residual term.

The α-stable distribution is highly intractable for inference because of the lack of a closed form density function in the general case. However, it is well-established that the α-stable distribution admits a Poisson series representation (PSR) in which the terms of the series are a function of the arrival times of a unit rate Poisson process. In our previous work, we have shown how to carry out inference for regression models using this series representation, which leads to a very convenient conditionally Gaussian framework, amenable to tractable Gaussian inference procedures. The PSR has to be truncated to a finite number of terms for practical purposes. The residual error terms have been approximated in our previous work by a Gaussian distribution, and we have recently shown that this approximation can be justified through a Central Limit Theorem (CLT). In this paper we present a new and exact characterisation of the first and second moments of the residual series over finite time intervals for the unit rate Poisson process, correcting a previous version that was only true in the infinite time limit. This enables us to test through simulation the rapid convergence of the residual terms to a Gaussian distribution of the Poisson series residual. We test this convergence using both Q-Q plots and the classical Kolmogorov-Smirnov test of Gaussianity.

In this paper we extend to the multidimensional case the modified Poisson series representation of linear stochastic processes driven by α-stable innovations. The latter has been recently introduced in the literature and it involves a Gaussian approximation of the residuals of the series, via the exact characterization of their moments. This allows for Bayesian techniques for parameter or state inference that would not be available otherwise, due to the lack of a closed-form likelihood function for the α-stable distribution. Simulation results are presented to validate the introduced extension and the quality of the approximation of the distribution. Finally, we show an example of generation from the process.

It is well known that the α-stable distribution, while having no closed form density function in the general case, admits a Poisson series representation (PSR) in which the terms of the series are a function of the arrival times of a unit rate Poisson process. In our previous work we have shown how to carry out inference for regression models using this series representation, which leads to a very convenient conditionally Gaussian framework, amenable to straightforward Gaussian inference procedures. The PSR has to be truncated to a finite number of terms for practical purposes. The residual terms have been approximated in our previous work by a Gaussian distribution with fully characterised moments. In this paper we present a new Central Limit Theorem (CLT) for the residual terms which serves to justify our previous approximation of the residual as Gaussian. Furthermore, we provide an analysis of the asymptotic convergence rate expressed in the CLT.

The α-stable distribution is very useful for modelling data with extreme values and skewed behaviour. The distribution is governed by two key parameters, tail thickness and skewness, in addition to scale and location. Inferring these parameters is difficult due to the lack of a closed form expression of the probability density. We develop a Bayesian method, based on the pseudo-marginal MCMC approach, that requires only unbiased estimates of the intractable likelihood. To compute these estimates we build an adaptive importance sampler for a latentvariable- representation of the α-stable density. This representation has previously been used in the literature for conditional MCMC sampling of the parameters, and we compare our method with this approach.

In this paper we develop an approach to Bayesian Monte Carlo inference for skewed α-stable distributions. Based on a series representation of the stable law in terms of infinite summations of random Poisson process arrival times, our framework leads to a simple representation in terms of conditionally Gaussian distributions for certain latent variables. Inference can therefore be carried out straightforwardly using techniques such as auxiliary variables versions of Markov chain Monte Carlo (MCMC) methods. The Poisson series representation (PSR) is further extended to practical application by introducing an approximation of the series residual terms based on exact moment calculations. Simulations illustrate the proposed framework applied to skewed α-stable simulated and real-world data, successfully estimating the distribution parameter values and being consistent with other (non-Bayesian) approaches. The methods are highly suitable for incorporation into hierarchical Bayesian models, and in this case the conditionally Gaussian structure of our model will lead to very efficient computations compared to other approaches.