High-dimensional Posterior Distribution Research Articles

Bayesian Estimation of Fault Slip Distribution for Slow Slip Events Based on an Efficient Hybrid Optimal Directional Gibbs Sampler and Its Application to the Guerrero 2006 Event

We propose an efficient Bayesian approach to infer a fault displacement from geodetic data in a slow slip event. Our physical model of the slip process reduces to a multiple linear regression subject to constraints. Assuming a Gaussian model for the geodetic data and considering a multivariate truncated normal prior distribution for the unknown fault slip, the resulting posterior distribution is also multivariate truncated normal. Regarding the posterior, we propose an algorithm based on Optimal Directional Gibbs that allows us to efficiently sample from the resulting high-dimensional posterior distribution of along dip and along strike movements of our fault grid division. A synthetic fault slip example illustrates the flexibility and accuracy of the proposed approach. The methodology is also applied to a real data set, for the 2006 Guerrero, Mexico, Slow Slip Event, where the objective is to recover the fault slip on a known interface that produces displacements observed at ground geodetic stations. As a by-product of our approach, we are able to estimate moment magnitude for the 2006 Guerrero Event with uncertainty quantification.

Identification of a monitoring nonlinear oil damper using particle filtering approach

Oil dampers have been used in recent years for passive structural control and shock mitigation in dynamic structural systems. However, machining technical and economic problems limit its application for new and existing buildings due to the necessary precision machining of incorporating shaft bearings and pressure sealings. A kind of nonlinear oil damper that is quite different from traditional oil damper is investigated and applied to a steel building. Vibration monitoring system was instrumented on the dampers to explore their actual performance and effectivity under strong earthquakes. Based on monitoring response of the nonlinear dampers under various excitations, Bayesian model selection is employed to analyze the most probable model class which can capture main dynamic characteristics of the nonlinear oil dampers and can also be used for predicting future response as well as reliability. Then, a particle filtering approach is proposed to identify the nonlinear model of the damper and quantify the model uncertainty. The developed particle filter is capable of re-parameterizing joint posterior distribution of states and parameters of the nonlinear oil damper without augmented state estimation, which combined with Markov chain Monte Carlo algorithm so as to be able to sample high-dimensional posterior distribution. The identified models and posterior distributions of parameters show that the developed particle filter approach can be appropriately used for nonlinear parameter identification without stuck to special particles. Furthermore, the dynamic properties of the nonlinear oil damper with respect to various excitations involving different spectral characteristics are discussed.

On polynomial-time computation of high-dimensional posterior measures by Langevin-type algorithms

The problem of generating random samples of high-dimensional posterior distributions is considered. The main results consist of non-asymptotic computational guarantees for Langevin-type MCMC algorithms which scale polynomially in key quantities such as the dimension of the model, the desired precision level, and the number of available statistical measurements. As a direct consequence, it is shown that posterior mean vectors as well as optimisation based maximum a posteriori (MAP) estimates are computable in polynomial time, with high probability under the distribution of the data. These results are complemented by statistical guarantees for recovery of the ground truth parameter generating the data. Our results are derived in a general high-dimensional non-linear regression setting (with Gaussian process priors) where posterior measures are not necessarily log-concave, employing a set of local ‘geometric’ assumptions on the parameter space, and assuming that a good initialiser of the algorithm is available. The theory is applied to a representative non-linear example from PDEs involving a steady-state Schrödinger equation.

VB-DeepONet: A Bayesian operator learning framework for uncertainty quantification

Neural network based data-driven operator learning schemes have shown tremendous potential in computational mechanics. DeepONet is one such neural network architecture which has gained widespread appreciation owing to its excellent prediction capabilities. Having said that, being set in a deterministic framework exposes DeepONet architecture to the risk of overfitting, poor generalization and in its unaltered form, it is incapable of quantifying the uncertainties associated with its predictions. To address these challenges, we propose a novel Bayesian operator learning framework referred to as the Variational Bayes DeepONet (VB-DeepONet). VB-DeepONet is rooted in Bayesian statistics and hence, (a) is less prone to overfitting as compared to its deterministic counterpart, (b) has better generalization, and (c) yields predictive uncertainty which is instrumental when decision making is concerned. VB-DeepONet exploits variational inference and hence has the capacity to take into account high dimensional posterior distributions while keeping the associated computational cost reasonable. Different examples covering mechanics problems like diffusion reaction, gravity pendulum, advection diffusion have been considered to illustrate the performance of the proposed VB-DeepONet and comparisons have been drawn against DeepONet set in deterministic framework, Proper Orthogonal Decomposition based Gaussian Process and DenseED. The results obtained illustrate the efficacy of the proposed approach in solving uncertainty quantification problems.

Robust sampling for weak lensing and clustering analyses with the Dark Energy Survey

ABSTRACTRecent cosmological analyses rely on the ability to accurately sample from high-dimensional posterior distributions. A variety of algorithms have been applied in the field, but justification of the particular sampler choice and settings is often lacking. Here, we investigate three such samplers to motivate and validate the algorithm and settings used for the Dark Energy Survey (DES) analyses of the first 3 yr (Y3) of data from combined measurements of weak lensing and galaxy clustering. We employ the full DES Year 1 likelihood alongside a much faster approximate likelihood, which enables us to assess the outcomes from each sampler choice and demonstrate the robustness of our full results. We find that the ellipsoidal nested sampling algorithm multinest reports inconsistent estimates of the Bayesian evidence and somewhat narrower parameter credible intervals than the sliced nested sampling implemented in polychord. We compare the findings from multinest and polychord with parameter inference from the Metropolis–Hastings algorithm, finding good agreement. We determine that polychord provides a good balance of speed and robustness for posterior and evidence estimation, and recommend different settings for testing purposes and final chains for analyses with DES Y3 data. Our methodology can readily be reproduced to obtain suitable sampler settings for future surveys.

Prior normalization for certified likelihood-informed subspace detection of Bayesian inverse problems

Markov chain Monte Carlo (MCMC) methods form one of the algorithmic foundations of Bayesian inverse problems. The recent development of likelihood-informed subspace (LIS) methods offers a viable route to designing efficient MCMC methods for exploring high-dimensional posterior distributions via exploiting the intrinsic low-dimensional structure of the underlying inverse problem. However, existing LIS methods and the associated performance analysis often assume that the prior distribution is Gaussian. This assumption is limited for inverse problems aiming to promote sparsity in the parameter estimation, as heavy-tailed priors, e.g., Laplace distribution or the elastic net commonly used in Bayesian LASSO, are often needed in this case. To overcome this limitation, we consider a prior normalization technique that transforms any non-Gaussian (e.g. heavy-tailed) priors into standard Gaussian distributions, which makes it possible to implement LIS methods to accelerate MCMC sampling via such transformations. We also rigorously investigate the integration of such transformations with several MCMC methods for high-dimensional problems. Finally, we demonstrate various aspects of our theoretical claims on two nonlinear inverse problems.

Approximate Bayesian neural Doppler imaging

Aims. The non-uniform surface temperature distribution of rotating active stars is routinely mapped with the Doppler imaging technique. Inhomogeneities in the surface produce features in high-resolution spectroscopic observations that shift in wavelength because of the Doppler effect, depending on their position on the visible hemisphere. The inversion problem has been systematically solved using maximum a posteriori regularized methods assuming smoothness or maximum entropy. Our aim in this work is to solve the full Bayesian inference problem by providing access to the posterior distribution of the surface temperature in the star compatible with the observations. Methods. We use amortized neural posterior estimation to produce a model that approximates the high-dimensional posterior distribution for spectroscopic observations of selected spectral ranges sampled at arbitrary rotation phases. The posterior distribution is approximated with conditional normalizing flows, which are flexible, tractable, and easy-to-sample approximations to arbitrary distributions. When conditioned on the spectroscopic observations, these normalizing flows provide a very efficient way of obtaining samples from the posterior distribution. The conditioning on observations is achieved through the use of Transformer encoders, which can deal with arbitrary wavelength sampling and rotation phases. Results. Our model can produce thousands of posterior samples per second, each one accompanied by an estimation of the log-probability. Our exhaustive validation of the model for very high-signal-to-noise observations shows that it correctly approximates the posterior, albeit with some overestimation of the broadening. We apply the model to the moderately fast rotator II Peg, producing the first Bayesian map of its temperature inhomogenities. We conclude that conditional normalizing flows are a very promising tool for carrying out approximate Bayesian inference in more complex problems in stellar physics, such as constraining the magnetic properties using polarimetry.

Bayesian non-linear regression with spatial priors for noise reduction and error estimation in quantitative MRI with an application in T1 estimation

Purpose. To develop a method that can reduce and estimate uncertainty in quantitative MR parameter maps without the need for hand-tuning of any hyperparameters. Methods. We present an estimation method where uncertainties are reduced by incorporating information on spatial correlations between neighbouring voxels. The method is based on a Bayesian hierarchical non-linear regression model, where the parameters of interest are sampled, using Markov chain Monte Carlo (MCMC), from a high-dimensional posterior distribution with a spatial prior. The degree to which the prior affects the model is determined by an automatic hyperparameter search using an information criterion and is, therefore, free from manual user-dependent tuning. The samples obtained further provide a convenient means to obtain uncertainties in both voxels and regions. The developed method was evaluated on T 1 estimations based on the variable flip angle method. Results. The proposed method delivers noise-reduced T 1 parameter maps with associated error estimates by combining MCMC sampling, the widely applicable information criterion, and total variation-based denoising. The proposed method results in an overall decrease in estimation error when compared to conventional voxel-wise maximum likelihood estimation. However, this comes with an increased bias in some regions, predominately at tissue interfaces, as well as an increase in computational time. Conclusions. This study provides a method that generates more precise estimates compared to the conventional method, without incorporating user subjectivity, and with the added benefit of uncertainty estimation.

Weighted Bayesian bootstrap for scalable posterior distributions

AbstractWe introduce and develop a weighted Bayesian bootstrap (WBB) for machine learning and statistics. WBB provides uncertainty quantification by sampling from a high dimensional posterior distribution. WBB is computationally fast and scalable using only off‐the‐shelf optimization software. First‐order asymptotic analysis provides a theoretical justification under suitable regularity conditions on the statistical model. We illustrate the proposed methodology in regularized regression, trend filtering and deep learning and conclude with directions for future research.

Approximation and sampling of multivariate probability distributions in the tensor train decomposition

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, a methodology that has been exploited for many years for scalable, high-dimensional density function approximation in quantum physics and chemistry. We build upon recent developments of the cross approximation algorithms in linear algebra to construct a tensor train approximation to the target probability density function using a small number of function evaluations. For sufficiently smooth distributions, the storage required for accurate tensor train approximations is moderate, scaling linearly with dimension. In turn, the structure of the tensor train surrogate allows sampling by an efficient conditional distribution method since marginal distributions are computable with linear complexity in dimension. Expected values of non-smooth quantities of interest, with respect to the surrogate distribution, can be estimated using transformed independent uniformly-random seeds that provide Monte Carlo quadrature or transformed points from a quasi-Monte Carlo lattice to give more efficient quasi-Monte Carlo quadrature. Unbiased estimates may be calculated by correcting the transformed random seeds using a Metropolis–Hastings accept/reject step, while the quasi-Monte Carlo quadrature 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; thus, the integrated autocorrelation time may be made arbitrarily close to 1, implying that, asymptotic in sample size, the cost per effectively independent sample is one target density evaluation plus the cheap tensor train surrogate proposal that has linear cost with dimension. 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. In all computed examples, the importance weight-corrected quasi-Monte Carlo quadrature performs best and is more efficient than DRAM by orders of magnitude across a wide range of approximation accuracies and sample sizes. Indeed, all the methods developed here significantly outperform DRAM in all computed examples.

Localization for MCMC: sampling high-dimensional posterior distributions with local structure

We investigate how ideas from covariance localization in numerical weather prediction can be used in Markov chain Monte Carlo (MCMC) sampling of high-dimensional posterior distributions arising in Bayesian inverse problems. To localize an inverse problem is to enforce an anticipated “local” structure by (i) neglecting small off-diagonal elements of the prior precision and covariance matrices; and (ii) restricting the influence of observations to their neighborhood. For linear problems we can specify the conditions under which posterior moments of the localized problem are close to those of the original problem. We explain physical interpretations of our assumptions about local structure and discuss the notion of high dimensionality in local problems, which is different from the usual notion of high dimensionality in function space MCMC. The Gibbs sampler is a natural choice of MCMC algorithm for localized inverse problems and we demonstrate that its convergence rate is independent of dimension for localized linear problems. Nonlinear problems can also be tackled efficiently by localization and, as a simple illustration of these ideas, we present a localized Metropolis-within-Gibbs sampler. Several linear and nonlinear numerical examples illustrate localization in the context of MCMC samplers for inverse problems.

The Recycling Gibbs sampler for efficient learning

Monte Carlo methods are essential tools for Bayesian inference. Gibbs sampling is a well-known Markov chain Monte Carlo (MCMC) algorithm, extensively used in signal processing, machine learning, and statistics, employed to draw samples from complicated high-dimensional posterior distributions. The key point for the successful application of the Gibbs sampler is the ability to draw efficiently samples from the full-conditional probability density functions. Since in the general case this is not possible, in order to speed up the convergence of the chain, it is required to generate auxiliary samples whose information is eventually disregarded. In this work, we show that these auxiliary samples can be recycled within the Gibbs estimators, improving their efficiency with no extra cost. This novel scheme arises naturally after pointing out the relationship between the standard Gibbs sampler and the chain rule used for sampling purposes. Numerical simulations involving simple and real inference problems confirm the excellent performance of the proposed scheme in terms of accuracy and computational efficiency. In particular we give empirical evidence of performance in a toy example, inference of Gaussian processes hyperparameters, and learning dependence graphs through regression.

Automated Sensitivity Analysis for Bayesian Inference via Markov Chain Monte Carlo: Applications to Gibbs Sampling

Bayesian inference relies heavily on numerical Markov chain Monte carlo (MCMC) methods for the estimation of the typically intractable high-dimensional posterior distributions and requires specific inputs. In this paper we introduce a new general and efficient numerical approach to address important robustness concerns of MCMC analysis with respect to prior input assumptions, a major obstacle to wider acceptance of Bayesian inference, including MCMC algorithm performance (convergence) reflected in the dependence on the chain starting values. The approach builds on recent developments in sensitivity analysis of high-dimensional numerical integrals for classical simulation methods using automatic numerical differentiation methods to compute first order derivatives of algorithmic output with respect to all inputs. We introduce a range of new robustness measures based on Jacobian matrices of MCMC output w.r.t. to the two sets of input parameters, prior parameters and chain starting values, to enable researchers to routinely undertake a comprehensive sensitivity analysis of their MCMC results. The methods are implemented for a range of Gibbs samplers and illustrated using both simulated and real data examples.

Fast Gibbs sampling for high-dimensional Bayesian inversion

Solving ill-posed inverse problems by Bayesian inference has recently attracted considerable attention. Compared to deterministic approaches, the probabilistic representation of the solution by the posterior distribution can be exploited to explore and quantify its uncertainties. In applications where the inverse solution is subject to further analysis procedures can be a significant advantage. Alongside theoretical progress, various new computational techniques allow us to sample very high dimensional posterior distributions: in (Lucka 2012 Inverse Problems 28 125012), and a Markov chain Monte Carlo posterior sampler was developed for linear inverse problems with -type priors. In this article, we extend this single component (SC) Gibbs-type sampler to a wide range of priors used in Bayesian inversion, such as general priors with additional hard constraints. In addition, a fast computation of the conditional, SC densities in an explicit, parameterized form, a fast, robust and exact sampling from these one-dimensional densities is key to obtain an efficient algorithm. We demonstrate that a generalization of slice sampling can utilize their specific structure for this task and illustrate the performance of the resulting slice-within-Gibbs samplers by different computed examples. These new samplers allow us to perform sample-based Bayesian inference in high-dimensional scenarios with certain priors for the first time, including the inversion of computed tomography data with the popular isotropic total variation prior.

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 adaptive SMC methods. The theoretical framework developed here will cover, under assumptions, several commonly used SMC algorithms [Chopin, Biometrika 89 (2002) 539-551; Jasra et al., Scand. J. Stat. 38 (2011) 1-22; Schafer and Chopin, Stat. Comput. 23 (2013) 163- 184]. 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 contexts [Jasra et al. (2011); Schafer and Chopin (2013)]. We establish that for a general class of adaptive SMC algorithms [Chopin (2002)], the asymptotic variance of the estimators from the adaptive SMC method is identical to a limiting 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 highdimensional parameters of the proposal kernels is critical for the efficiency of the algorithm.

Bayesian analysis of circular distributions based on non-negative trigonometric sums

ABSTRACTFernández-Durán [Circular distributions based on nonnegative trigonometric sums. Biometrics. 2004;60:499–503] developed a new family of circular distributions based on non-negative trigonometric sums that is suitable for modelling data sets that present skewness and/or multimodality. In this paper, a Bayesian approach to deriving estimates of the unknown parameters of this family of distributions is presented. Because the parameter space is the surface of a hypersphere and the dimension of the hypersphere is an unknown parameter of the distribution, the Bayesian inference must be based on transdimensional Markov Chain Monte Carlo (MCMC) algorithms to obtain samples from the high-dimensional posterior distribution. The MCMC algorithm explores the parameter space by moving along great circles on the surface of the hypersphere. The methodology is illustrated with real and simulated data sets.

The VIMOS Public Extragalactic Redshift Survey

Aims. Using the VIMOS Public Extragalactic Redshift Survey (VIPERS) we aim to jointly estimate the key parameters that describe the galaxy density field and its spatial correlations in redshift space. Methods. We use the Bayesian formalism to jointly reconstruct the redshift-space galaxy density field, power spectrum, galaxy bias and galaxy luminosity function given the observations and survey selection function. The high-dimensional posterior distribution is explored using the Wiener filter within a Gibbs sampler. We validate the analysis using simulated catalogues and apply it to VIPERS data taking into consideration the inhomogeneous selection function. Results. We present joint constraints on the anisotropic power spectrum as well as the bias and number density of red and blue galaxy classes in luminosity and redshift bins as well as the measurement covariances of these quantities. We find that the inferred galaxy bias and number density parameters are strongly correlated although these are only weakly correlated with the galaxy power spectrum. The power spectrum and redshift-space distortion parameters are in agreement with previous VIPERS results with the value of the growth rate $f\sigma_8 = 0.38$ with 18% uncertainty at redshift 0.7.

Sequential Monte Carlo for Maximum Weight Subgraphs with Application to Solving Image Jigsaw Puzzles.

We consider a problem of finding maximum weight subgraphs (MWS) that satisfy hard constraints in a weighted graph. The constraints specify the graph nodes that must belong to the solution as well as mutual exclusions of graph nodes, i.e., pairs of nodes that cannot belong to the same solution. Our main contribution is a novel inference approach for solving this problem in a sequential monte carlo (SMC) sampling framework. Usually in an SMC framework there is a natural ordering of the states of the samples. The order typically depends on observations about the states or on the annealing setup used. In many applications (e.g., image jigsaw puzzle problems), all observations (e.g., puzzle pieces) are given at once and it is hard to define a natural ordering. Therefore, we relax the assumption of having ordered observations about states and propose a novel SMC algorithm for obtaining maximum a posteriori estimate of a high-dimensional posterior distribution. This is achieved by exploring different orders of states and selecting the most informative permutations in each step of the sampling. Our experimental results demonstrate that the proposed inference framework significantly outperforms loopy belief propagation in solving the image jigsaw puzzle problem. In particular, our inference quadruples the accuracy of the puzzle assembly compared to that of loopy belief propagation.

Computational statistics using the Bayesian Inference Engine

This paper introduces the Bayesian Inference Engine (BIE), a general parallel, optimised software package for parameter inference and model selection. This package is motivated by the analysis needs of modern astronomical surveys and the need to organise and reuse expensive derived data. The BIE is the first platform for computational statistics designed explicitly to enable Bayesian update and model comparison for astronomical problems. Bayesian update is based on the representation of high-dimensional posterior distributions using metric-ball-tree based kernel density estimation. Among its algorithmic offerings, the BIE emphasises hybrid tempered MCMC schemes that robustly sample multimodal posterior distributions in high-dimensional parameter spaces. Moreover, the BIE is implements a full persistence or serialisation system that stores the full byte-level image of the running inference and previously characterised posterior distributions for later use. Two new algorithms to compute the marginal likelihood from the posterior distribution, developed for and implemented in the BIE, enable model comparison for complex models and data sets. Finally, the BIE was designed to be a collaborative platform for applying Bayesian methodology to astronomy. It includes an extensible object-oriented and easily extended framework that implements every aspect of the Bayesian inference. By providing a variety of statistical algorithms for all phases of the inference problem, a scientist may explore a variety of approaches with a single model and data implementation. Additional technical details and download details are available from http://www.astro.umass.edu/bie. The BIE is distributed under the GNU GPL.

A New Bayesian Method Incorporating With Local Correlation for IBM Estimation

A lot of efforts have been made in the Ideal Binary Mask (IBM) estimation via statistical learning methods. The Bayesian method is a common one. However, one drawback is that the mask is estimated for each time-frequency (T-F) unit independently. The correlation between units has not been fully taken into account. In this paper, we attempt to consider the local correlation information from two aspects to improve the performance. On one hand, a T-F segmentation based potential function is proposed to depict the local correlation between the mask labels of adjacent units directly. It is derived from a demonstrated assumption that units which belong to one segment are mainly dominated by one source. On the other hand, a local noise level tracking stage is incorporated. The local level is obtained by averaging among several adjacent units and can be considered as an approach to true noise energy. It is used as the intermediary auxiliary variable to indicate the correlation. While some secondary factors are omitted, the high dimensional posterior distribution is simulated by a Markov Chain Monte Carlo (MCMC) method. In iterations, the correlation is fully considered to compute the acceptance ratio. The estimate of IBM is obtained by the expectation. Our system is evaluated and compared with previous Bayesian system, and it yields substantially better performance in terms of HIT-FA rates and SNR gain.