Search

In this chapter, we address the challenge of exploring the posterior distributions of Bayesian inverse problems with computationally intensive forward models. We consider various multivariate proposal distributions, and compare them with single-site Metropolis updates. We show how fast, approximate models can be leveraged to improve the MCMC sampling efficiency., Comment: To appear in the Handbook of Markov Chain Monte Carlo (2nd edition)

We develop a novel Markov chain Monte Carlo (MCMC) method that exploits a hierarchy of models of increasing complexity to efficiently generate samples from an unnormalized target distribution. Broadly, the method rewrites the Multilevel MCMC approach of Dodwell et al. (2015) in terms of the Delayed Acceptance (DA) MCMC of Christen & Fox (2005). In particular, DA is extended to use a hierarchy of models of arbitrary depth, and allow subchains of arbitrary length. We show that the algorithm satisfies detailed balance, hence is ergodic for the target distribution. Furthermore, multilevel variance reduction is derived that exploits the multiple levels and subchains, and an adaptive multilevel correction to coarse-level biases is developed. Three numerical examples of Bayesian inverse problems are presented that demonstrate the advantages of these novel methods. The software and examples are available in PyMC3., Comment: 29 pages, 12 figures

The standard models of sequence evolution on a tree determine probabilities for every character or site pattern. A flattening is an arrangement of these probabilities into a matrix, with rows corresponding to all possible site patterns for one set $A$ of taxa and columns corresponding to all site patterns for another set $B$ of taxa. Flattenings have been used to prove difficult results relating to phylogenetic invariants and consistency and also form the basis of several methods of phylogenetic inference. We prove that the rank of the flattening equals $r^{\ell_T(A|B)}$, where $r$ is the number of states and $\ell_T(A|B)$ is the parsimony length of the binary character separating $A$ and $B$. This result corrects an earlier published formula and opens up new applications for old parsimony theorems. Since completing this work, we have learnt that an equivalent result has been proved much earlier by Casanellas and Fern\'andez-S\'anchez, using a different proof strategy.

We address the inverse Frobenius--Perron problem: given a prescribed target distribution $\rho$, find a deterministic map $M$ such that iterations of $M$ tend to $\rho$ in distribution. We show that all solutions may be written in terms of a factorization that combines the forward and inverse Rosenblatt transformations with a uniform map, that is, a map under which the uniform distribution on the $d$-dimensional hypercube as invariant. Indeed, every solution is equivalent to the choice of a uniform map. We motivate this factorization via $1$-dimensional examples, and then use the factorization to present solutions in $1$ and $2$ dimensions induced by a range of uniform maps., Comment: Submitted to Entropy. v2: a few typos fixed

Uncertainty Quantification through Markov Chain Monte Carlo (MCMC) can be prohibitively expensive for target probability densities with expensive likelihood functions, for instance when the evaluation it involves solving a Partial Differential Equation (PDE), as is the case in a wide range of engineering applications. Multilevel Delayed Acceptance (MLDA) with an Adaptive Error Model (AEM) is a novel approach, which alleviates this problem by exploiting a hierarchy of models, with increasing complexity and cost, and correcting the inexpensive models on-the-fly. The method has been integrated within the open-source probabilistic programming package PyMC3 and is available in the latest development version. In this paper, the algorithm is presented along with an illustrative example., Comment: 8 pages, 4 figures, accepted for Machine Learning for Engineering Modeling, Simulation, and Design Workshop at Neural Information Processing Systems 2020

Bayesian modelling and computational inference by Markov chain Monte Carlo (MCMC) is a principled framework for large-scale uncertainty quantification, though is limited in practice by computational cost when implemented in the simplest form that requires simulating an accurate computer model at each iteration of the MCMC. The delayed acceptance Metropolis--Hastings MCMC leverages a reduced model for the forward map to lower the compute cost per iteration, though necessarily reduces statistical efficiency that can, without care, lead to no reduction in the computational cost of computing estimates to a desired accuracy. Randomizing the reduced model for the forward map can dramatically improve computational efficiency, by maintaining the low cost per iteration but also avoiding appreciable loss of statistical efficiency. Randomized maps are constructed by a posteriori adaptive tuning of a randomized and locally-corrected deterministic reduced model. Equivalently, the approximated posterior distribution may be viewed as induced by a modified likelihood function for use with the reduced map, with parameters tuned to optimize the quality of the approximation to the correct posterior distribution. Conditions for adaptive MCMC algorithms allow practical approximations and algorithms that have guaranteed ergodicity for the target distribution. Good statistical and computational efficiencies are demonstrated in examples of calibration of large-scale numerical models of geothermal reservoirs and electrical capacitance tomography., Comment: This is a pre-print of an article submitted to GEM - International Journal on Geomathematics. arXiv admin note: text overlap with arXiv:1809.03176

We describe a new and computationally efficient Bayesian methodology for inferring species trees and demographics from unlinked binary markers. Likelihood calculations are carried out using diffusion models of allele frequency dynamics combined with a new algorithm for numerically computing likelihoods of quantitative traits. The diffusion approach allows for analysis of datasets containing hundreds or thousands of individuals. The method, which we call \snapper, has been implemented as part of the Beast2 package. We introduce the models, the efficient algorithms, and report performance of \snapper on simulated data sets and on SNP data from rattlesnakes and freshwater turtles.

General multivariate distributions are notoriously expensive to sample from, particularly the high-dimensional posterior distributions in PDE-constrained inverse problems. This paper develops a sampler for arbitrary continuous multivariate distributions that is based on low-rank surrogates in the tensor-train format. We construct a tensor-train approximation to the target probability density function using the cross interpolation, which requires a small number of function evaluations. For sufficiently smooth distributions the storage required for the TT approximation is moderate, scaling linearly with dimension. The structure of the tensor-train surrogate allows efficient sampling by the conditional distribution method. Unbiased estimates may be calculated by correcting the transformed random seeds using a Metropolis--Hastings accept/reject step. Moreover, one can use a more efficient quasi-Monte Carlo quadrature that may be corrected either by a control-variate strategy, or by importance weighting. We show that the error in the tensor-train approximation propagates linearly into the Metropolis--Hastings rejection rate and the integrated autocorrelation time of the resulting Markov chain. These methods are demonstrated in three computed examples: fitting failure time of shock absorbers; a PDE-constrained inverse diffusion problem; and sampling from the Rosenbrock distribution. The delayed rejection adaptive Metropolis (DRAM) algorithm is used as a benchmark. We find that the importance-weight corrected quasi-Monte Carlo quadrature performs best in all computed examples, and is orders-of-magnitude more efficient than DRAM across a wide range of approximation accuracies and sample sizes. Indeed, all the methods developed here significantly outperform DRAM in all computed examples., Comment: 32 pages

Sample-based Bayesian inference provides a route to uncertainty quantification in the geosciences, and inverse problems in general, though is very computationally demanding in the naive form that requires simulating an accurate computer model at each iteration. We present a new approach that constructs a stochastic correction to the error induced by a reduced model, with the correction improving as the algorithm proceeds. This enables sampling from the correct target distribution at reduced computational cost per iteration, as in existing delayed-acceptance schemes, while avoiding appreciable loss of statistical efficiency that necessarily occurs when using a reduced model. Use of the stochastic correction significantly reduces the computational cost of estimating quantities of interest within desired uncertainty bounds. In contrast, existing schemes that use a reduced model directly as a surrogate do not actually improve computational efficiency in our target applications. We build on recent simplified conditions for adaptive Markov chain Monte Carlo algorithms to give practical approximation schemes and algorithms with guaranteed convergence. The efficacy of this new approach is demonstrated in two computational examples, including calibration of a large-scale numerical model of a real geothermal reservoir, that show good computational and statistical efficiencies on both synthetic and measured data sets., Comment: Submitted to IJNME

Trajectory reconstruction is the process of inferring the path of a moving object between successive observations. In this paper, we propose a smoothing spline -- which we name the V-spline -- that incorporates position and velocity information and a penalty term that controls acceleration. We introduce a particular adaptive V-spline designed to control the impact of irregularly sampled observations and noisy velocity measurements. A cross-validation scheme for estimating the V-spline parameters is given and we detail the performance of the V-spline on four particularly challenging test datasets. Finally, an application of the V-spline to vehicle trajectory reconstruction in two dimensions is given, in which the penalty term is allowed to further depend on known operational characteristics of the vehicle., Comment: 25 pages, submitted

We consider posterior sampling in the very common Bayesian hierarchical model in which observed data depends on high-dimensional latent variables that, in turn, depend on relatively few hyperparameters. When the full conditional over the latent variables has a known form, the marginal posterior distribution over hyperparameters is accessible and can be sampled using a Markov chain Monte Carlo (MCMC) method on a low-dimensional parameter space. This may improve computational efficiency over standard Gibbs sampling since computation is not over the high-dimensional space of latent variables and correlations between hyperparameters and latent variables become irrelevant. When the marginal posterior over hyperparameters depends on a fixed-dimensional sufficient statistic, precomputation of the sufficient statistic renders the cost of the low-dimensional MCMC independent of data size. Then, when the hyperparameters are the primary variables of interest, inference may be performed in big-data settings at modest cost. Moreover, since the form of the full conditional for the latent variables does not depend on the form of the hyperprior distribution, the method imposes no restriction on the hyperprior, unlike Gibbs sampling that typically requires conjugate distributions. We demonstrate these efficiency gains in four computed examples., Comment: 24 pages, 3 figures

We develop a finite-dimensional approximation of the Frobenius-Perron operator using the finite volume method applied to the continuity equation for the evolution of probability. A Courant-Friedrichs-Lewy condition ensures that the approximation satisfies the Markov property, while existing convergence theory for the finite volume method guarantees convergence of the discrete operator to the continuous operator as mesh size tends to zero. Properties of the approximation are demonstrated in a computed example of sequential inference for the state of a low-dimensional mechanical system when observations give rise to multi-modal distributions., Comment: 15 pages, 2 figures

Proposals for Metropolis-Hastings MCMC derived by discretizing Langevin diffusion or Hamiltonian dynamics are examples of stochastic autoregressive proposals that form a natural wider class of proposals with equivalent computability. We analyze Metropolis-Hastings MCMC with stochastic autoregressive proposals applied to target distributions that are absolutely continuous with respect to some Gaussian distribution to derive expressions for expected acceptance probability and expected jump size, as well as measures of computational cost, in the limit of high dimension. Thus, we are able to unify existing analyzes for these classes of proposals, and to extend the theoretical results that provide useful guidelines for tuning the proposals for optimal computational efficiency. For the simplified Langevin algorithm we find that it is optimal to take at least three steps of the proposal before the Metropolis-Hastings accept-reject step, and for Hamiltonian/hybrid Monte Carlo we provide new guidelines for the optimal number of integration steps and criteria for choosing the optimal mass matrix., Comment: 44 pages, 1 figure. arXiv admin note: text overlap with arXiv:1605.05441

We analyse computational efficiency of Metropolis-Hastings algorithms with stochastic AR(1) process proposals. These proposals include, as a subclass, discretized Langevin diffusion (e.g. MALA) and discretized Hamiltonian dynamics (e.g. HMC). We derive expressions for the expected acceptance rate and expected jump size for MCMC methods with general stochastic AR(1) process proposals for the case where the target distribution is absolutely continuous with respect to a Gaussian and the covariance of the Gaussian is allowed to have off-diagonal terms. This allows us to extend what is known about several MCMC methods as well as determining the efficiency of new MCMC methods of this type. In the special case of Hybrid Monte Carlo, we can determine the optimal integration time and the effect of the choice of mass matrix. By including the effect of Metropolis-Hastings we also extend results by Fox and Parker, who used matrix splitting techniques to analyse the performance and improve efficiency of stochastic AR(1) processes for sampling from Gaussian distributions., Comment: 58 pages. Electronics Technical Reports No. 2016-1, ISSN 1172-496X, May 2016. arXiv admin note: substantial text overlap with arXiv:1501.03150

We present an efficient and flexible method for computing likelihoods of phenotypic traits on a phylogeny. The method does not resort to Monte-Carlo computation but instead blends Felsenstein's discrete character pruning algorithm with methods for numerical quadrature. It is not limited to Gaussian models and adapts readily to model uncertainty in the observed trait values. We demonstrate the framework by developing efficient algorithms for likelihood calculation and ancestral state reconstruction under Wright's threshold model, applying our methods to a dataset of trait data for extrafloral nectaries (EFNs) across a phylogeny of 839 Labales species., Comment: 23 pages, 3 figures

We solve the inverse problem of deblurring a pixelized image of Jupiter using regularized deconvolution and by sample-based Bayesian inference. By efficiently sampling the marginal posterior distribution for hyperparameters, then the full conditional for the deblurred image, we find that we can evaluate the posterior mean faster than regularized inversion, when selection of the regularizing parameter is considered. To our knowledge, this is the first demonstration of sampling and inference that takes less compute time than regularized inversion in an inverse problems. Comparison to random-walk Metropolis-Hastings and block Gibbs MCMC shows that marginal then conditional sampling also outperforms these more common sampling algorithms, having better scaling with problem size. When problem-specific computations are feasible the asymptotic cost of an independent sample is one linear solve, implying that sample-based Bayesian inference may be performed directly over function spaces, when that limit exists.

Standard Gibbs sampling applied to a multivariate normal distribution with a specified precision matrix is equivalent in fundamental ways to the Gauss-Seidel iterative solution of linear equations in the precision matrix. Specifically, the iteration operators, the conditions under which convergence occurs, and geometric convergence factors (and rates) are identical. These results hold for arbitrary matrix splittings from classical iterative methods in numerical linear algebra giving easy access to mature results in that field, including existing convergence results for antithetic-variable Gibbs sampling, REGS sampling, and generalizations. Hence, efficient deterministic stationary relaxation schemes lead to efficient generalizations of Gibbs sampling. The technique of polynomial acceleration that significantly improves the convergence rate of an iterative solver derived from a \emph{symmetric} matrix splitting may be applied to accelerate the equivalent generalized Gibbs sampler. Identicality of error polynomials guarantees convergence of the inhomogeneous Markov chain, while equality of convergence factors ensures that the optimal solver leads to the optimal sampler. Numerical examples are presented, including a Chebyshev accelerated SSOR Gibbs sampler applied to a stylized demonstration of low-level Bayesian image reconstruction in a large 3-dimensional linear inverse problem., Comment: 33 pages, 4 figures

We analyse computational efficiency of Metropolis-Hastings algorithms with AR(1) process proposals. These proposals include, as a subclass, discretized Langevin diffusion (e.g. MALA) and discretized Hamiltonian dynamics (e.g. HMC). By including the effect of Metropolis-Hastings we extend earlier work by Fox and Parker, who used matrix splitting techniques to analyse the performance and improve efficiency of AR(1) processes for targeting Gaussian distributions. Our research enables analysis of MCMC methods that draw samples from non-Gaussian target distributions by using AR(1) process proposals in Metropolis-Hastings algorithms, by analysing the matrix splitting of the precision matrix for a local Gaussian approximation of the non-Gaussian target., Comment: 34 pages

We consider Metropolis Hastings MCMC in cases where the log of the ratio of target distributions is replaced by an estimator. The estimator is based on m samples from an independent online Monte Carlo simulation. Under some conditions on the distribution of the estimator the process resembles Metropolis Hastings MCMC with a randomized transition kernel. When this is the case there is a correction to the estimated acceptance probability which ensures that the target distribution remains the equilibrium distribution. The simplest versions of the Penalty Method of Ceperley and Dewing (1999), the Universal Algorithm of Ball et al. (2003) and the Single Variable Exchange algorithm of Murray et al. (2006) are special cases. In many applications of interest the correction terms cannot be computed. We consider approximate versions of the algorithms. We show that on average O(m) of the samples realized by a simulation approximating a randomized chain of length n are exactly the same as those of a coupled (exact) randomized chain. Approximation biases Monte Carlo estimates with terms O(1/m) or smaller. This should be compared to the Monte Carlo error which is O(1/sqrt(n))., Comment: 20 pages, 5 graphs in 3 figures

Generalized Gibbs kernels are those that may take any direction not necessarily bounded to each axis along the parameters of the objective function. We study how to optimally choose such directions in a Directional, random scan, Gibbs sampler setting. The optimal direction is chosen by minimizing to the mutual information (Kullback-Leibler divergence) of two steps of the MCMC for a truncated Normal objective function. The result is generalized to be used when a Multivariate Normal (local) approximation is available for the objective function. Three Gibbs direction distributions are tested in highly skewed non-normal objective functions.