Exact Posterior Distribution Research Articles

Generating Independent Replicates Directly from the Posterior Distribution for a Class of Spatial Hierarchical Models

Markov chain Monte Carlo (MCMC) allows one to generate dependent replicates from a posterior distribution for effectively any Bayesian hierarchical model. However, MCMC can produce a significant computational burden. This motivates us to consider finding expressions of the posterior distribution that are computationally straightforward to obtain independent replicates from directly. We focus on a broad class of Bayesian hierarchical models for spatially dependent data, which are often modeled via a latent Gaussian process (LGP). First, we derive a new class of distributions referred to as the generalized conjugate multivariate (GCM) distribution. The GCM distribution’s theoretical development follows that of the conjugate multivariate (CM) distribution with two main differences: the GCM allows for latent Gaussian process assumptions, and the GCM explicitly accounts for hyperparameters through marginalization. The development of GCM is needed to obtain independent replicates directly from the exact posterior distribution, which has an efficient regression form. Hence, we refer to our method as Exact Posterior Regression (EPR). Simulation studies with weakly stationary spatial processes and spatial basis function expansions are provided. We provide an analysis of poverty incidence from the U.S. Census Bureau, and an analysis of high-dimensional remote sensing data. Supplementary materials for this article are available online.

An efficient method of predicting S-wave velocity using sparse Gaussian process regression for a tight sandstone reservoir

The shear wave (S-wave) velocity plays a crucial role in interpreting the lithology in seismic data, identifying fluids and predicting reservoirs. However, S-wave velocity is often unavailable due to the high cost of measurement and technical constraints. Conventional methods exhibit limitations that potentially impact the accuracy or efficiency on predicting S-wave velocity. Moreover, these methods always ignore the uncertainty quantification associated with the predicted results. This paper proposes a sparse Gaussian process regression (SGPR) method to predict the S-wave velocity in tight sandstone reservoirs. SGPR is a highly efficient regression technique that is based on the Gaussian process regression (GPR) method. In the SGPR method, inducing inputs are introduced to approximate the kernel matrix to decrease the computational complexity. A sparse set of inducing inputs and kernel hyperparameters are optimized through minimizing the Kullback-Leibler (KL) divergence between the exact posterior distribution and the approximate one. In this study, we select several types of logging data, which include porosity, water saturation, shale content, lithology and P-wave velocity, as the inputs for the SGPR method to predict S-wave velocity. To validate its effectiveness, we use the SGPR method to predict S-wave velocity in tight sandstone and compare the results with those from the GPR method, the bidirectional long short-term memory (BiLSTM) method and the Xu-White model. Additionally, we conduct cross-validation to demonstrate the robustness of the SGPR method. Our findings indicate that the SGPR method presents better performance and significant advantages about the accuracy and efficiency. Moreover, the SGPR method offers uncertainty quantification for the predicted S-wave velocity.

Exact Posterior distribution of risk ratio in the Kumaraswamy–Binomial model

Exact Posterior distribution of risk ratio in the Kumaraswamy–Binomial model

Wireless Beacon Enabled Hybrid Sparse Channel Estimation for RIS-Aided mmWave Communications

The inability of reconfigurable intelligent surfaces (RISs) to signal processing has become a critical bottleneck in terms of the precise capture of cascaded channels. There exists a pair of critical issues: i) the high-dimensional cascaded channel to be estimated, which always entails substantial pilot overhead typically proportional to RIS elements, and ii) the high possibilities of false alarms and unfaithful angle detection during cascaded channel estimates, referred to as the angular distortion, particularly in the case of conventional fully-passive RIS systems. To circumvent these issues, we investigate in this paper a multi-user multiple-input multiple-output (MIMO) system assisted by a wireless beacon (WB) enabled RIS over millimeter wave (mmWave) frequency. Specifically, WB is a functionally independent hardware that attaches parasitically to the RIS using only one radio frequency chain for broadcasting pilots, allowing for the accurate acquisition of partial cascaded channel information, e.g., some of angles of arrival/departure (AoA/AoDs) and complex gains. This efficiently eliminates the angular distortion effect and thus facilitates high precision recovery of sparse cascaded mmWave channels. Then, a hybrid structured sparsity (HSS) model is presented to capture hybrid sparse priors and approximate exact posterior distributions associated with RIS-aided twin-hop channels. We conceive of channel estimation as a compressive sensing problem in contrast to its conventional copy due to the presence of an unknown sensing matrix. To tackle this problem, an expectation-maximization (EM)-based algorithm, i.e., HSS-EM, is developed in order to fully employ the sparse priors encapsulated by the proposed HSS model for a robust and accurate recovery of the cascaded channels. Numerical results validate our analytical claims and demonstrate the significant improvement in terms of the normalized mean square error (NMSE) performance over conventional designs.

Obesity and labour market outcomes in Italy: a dynamic panel data evidence with correlated random effects.

This paper investigates the effects of obesity, socio-economic variables, and individual-specific factors on work productivity across Italian regions. A dynamic panel data with correlated random effects is used to jointly deal with incidental parameters, endogeneity issues, and functional forms of misspecification. Methodologically, a hierarchical semiparametric Bayesian approach is involved in shrinking high dimensional model classes, and then obtaining a subset of potential predictors affecting outcomes. Monte Carlo designs are addressed to construct exact posterior distributions and then perform accurate forecasts. Cross-sectional Heterogeneity is modelled nonparametrically allowing for correlation between heterogeneous parameters and initial conditions as well as individual-specific regressors. Prevention policies and strategies to handle health and labour market prospects are also discussed.

Fast, Scalable Approximations to Posterior Distributions in Extended Latent Gaussian Models

ABSTRACT We define a novel class of additive models, called Extended Latent Gaussian Models, that allow for a wide range of response distributions and flexible relationships between the additive predictor and mean response. The new class covers a broad range of interesting models including multi-resolution spatial processes, partial likelihood-based survival models, and multivariate measurement error models. Because computation of the exact posterior distribution is infeasible, we develop a fast, scalable approximate Bayesian inference methodology for this class based on nested Gaussian, Laplace, and adaptive quadrature approximations. We prove that the error in these approximate posteriors is under standard conditions, and provide numerical evidence suggesting that our method runs faster and scales to larger datasets than methods based on Integrated Nested Laplace Approximations and Markov chain Monte Carlo, with comparable accuracy. We apply the new method to the mapping of malaria incidence rates in continuous space using aggregated data, mapping leukemia survival hazards using a Cox Proportional-Hazards model with a continuously-varying spatial process, and estimating the mass of the Milky Way Galaxy using noisy multivariate measurements of the positions and velocities of star clusters in its orbit. Supplementary materials for this article are available online.

Model Reduction of Linear Dynamical Systems via Balancing for Bayesian Inference

We consider the Bayesian approach to the linear Gaussian inference problem of inferring the initial condition of a linear dynamical system from noisy output measurements taken after the initial time. In practical applications, the large dimension of the dynamical system state poses a computational obstacle to computing the exact posterior distribution. Model reduction offers a variety of computational tools that seek to reduce this computational burden. In particular, balanced truncation is a system-theoretic approach to model reduction which obtains an efficient reduced-dimension dynamical system by projecting the system operators onto state directions which trade off the reachability and observability of state directions as expressed through the associated Gramians. We introduce Gramian definitions relevant to the inference setting and propose a balanced truncation approach based on these inference Gramians that yield a reduced dynamical system that can be used to cheaply approximate the posterior mean and covariance. Our definitions exploit natural connections between (i) the reachability Gramian and the prior covariance and (ii) the observability Gramian and the Fisher information. The resulting reduced model then inherits stability properties and error bounds from system theoretic considerations, and in some settings yields an optimal posterior covariance approximation. Numerical demonstrations on two benchmark problems in model reduction show that our method can yield near-optimal posterior covariance approximations with order-of-magnitude state dimension reduction.

Bayesian Inference on Multivariate Medians and Quantiles

In this paper, we consider Bayesian inference on a class of multivariate median and the multivariate quantile functionals of a joint distribution using a Dirichlet process prior. Since, unlike univariate quantiles, the exact posterior distribution of multivariate median and multivariate quantiles are not obtainable explicitly, we study these distributions asymptotically. We derive a Bernstein-von Mises theorem for the multivariate $\ell_1$-median with respect to general $\ell_p$-norm, which in particular shows that its posterior concentrates around its true value at $n^{-1/2}$-rate and its credible sets have asymptotically correct frequentist coverage. In particular, asymptotic normality results for the empirical multivariate median with general $\ell_p$-norm is also derived in the course of the proof which extends the results from the case $p=2$ in the literature to a general $p$. The technique involves approximating the posterior Dirichlet process by a Bayesian bootstrap process and deriving a conditional Donsker theorem. We also obtain analogous results for an affine equivariant version of the multivariate $\ell_1$-median based on an adaptive transformation and re-transformation technique. The results are extended to a joint distribution of multivariate quantiles. The accuracy of the asymptotic result is confirmed by a simulation study. We also use the results to obtain Bayesian credible regions for multivariate medians for Fisher's iris data, which consists of four features measured for each of three plant species.

Ultralarge-scale approximations and galaxy clustering: Debiasing constraints on cosmological parameters

ABSTRACT Upcoming galaxy surveys will allow us to probe the growth of the cosmic large-scale structure with improved sensitivity compared to current missions, and will also map larger areas of the sky. This means that in addition to the increased precision in observations, future surveys will also access the ultralarge-scale regime, where commonly neglected effects such as lensing, redshift-space distortions, and relativistic corrections become important for calculating correlation functions of galaxy positions. At the same time, several approximations usually made in these calculations such as the Limber approximation break down at those scales. The need to abandon these approximations and simplifying assumptions at large scales creates severe issues for parameter estimation methods. On the one hand, exact calculations of theoretical angular power spectra become computationally expensive, and the need to perform them thousands of times to reconstruct posterior probability distributions for cosmological parameters makes the approach unfeasible. On the other hand, neglecting relativistic effects and relying on approximations may significantly bias the estimates of cosmological parameters. In this work, we quantify this bias and investigate how an incomplete modelling of various effects on ultralarge scales could lead to false detections of new physics beyond the standard ΛCDM model. Furthermore, we propose a simple debiasing method that allows us to recover true cosmologies without running the full parameter estimation pipeline with exact theoretical calculations. This method can therefore provide a fast way of obtaining accurate values of cosmological parameters and estimates of exact posterior probability distributions from ultralarge-scale observations.

Emulation-accelerated Hamiltonian Monte Carlo algorithms for parameter estimation and uncertainty quantification in differential equation models

We propose to accelerate Hamiltonian and Lagrangian Monte Carlo algorithms by coupling them with Gaussian processes for emulation of the log unnormalised posterior distribution. We provide proofs of detailed balance with respect to the exact posterior distribution for these algorithms, and validate the correctness of the samplers’ implementation by Geweke consistency tests. We implement these algorithms in a delayed acceptance (DA) framework, and investigate whether the DA scheme can offer computational gains over the standard algorithms. A comparative evaluation study is carried out to assess the performance of the methods on a series of models described by differential equations, including a real-world application of a 1D fluid-dynamics model of the pulmonary blood circulation. The aim is to identify the algorithm which gives the best trade-off between accuracy and computational efficiency, to be used in nonlinear DE models, which are computationally onerous due to repeated numerical integrations in a Bayesian analysis. Results showed no advantage of the DA scheme over the standard algorithms with respect to several efficiency measures based on the effective sample size for most methods and DE models considered. These gradient-driven algorithms register a high acceptance rate, thus the number of expensive forward model evaluations is not significantly reduced by the first emulator-based stage of DA. Additionally, the Lagrangian Dynamical Monte Carlo and Riemann Manifold Hamiltonian Monte Carlo tended to register the highest efficiency (in terms of effective sample size normalised by the number of forward model evaluations), followed by the Hamiltonian Monte Carlo, and the No U-turn sampler tended to be the least efficient.

Nonnegative Bayesian nonparametric factor models with completely random measures

We present a Bayesian nonparametric Poisson factorization model for modeling dense network data with an unknown and potentially growing number of overlapping communities. The construction is based on completely random measures and allows the number of communities to either increase with the number of nodes at a specified logarithmic or polynomial rate, or be bounded. We develop asymptotics for the number and size of the communities of the network and derive a Markov chain Monte Carlo algorithm for targeting the exact posterior distribution for this model. The usefulness of the approach is illustrated on various real networks.

Bayesian Inference in Common Microeconometric Models With Massive Datasets by Double Marginalized Subsampling

Bayesian inference with a large dataset is computationally intensive, as Markov chain Monte Carlo simulation requires a complete scan of the dataset for each proposed parameter update. To reduce the number of data points evaluated at each iteration of posterior simulation, we develop a double marginalized subsampling method, which is applicable to a wide array of microeconometric models including Tobit, Probit, regressions with non-Gaussian errors, heteroscedasticity and stochastic volatility, hierarchical longitudinal models, time-varying-parameter regressions, Gaussian mixtures, etc. We also provide an extension to double pseudo-marginalized subsampling, which has more applications beyond conditionally conjugate models. With rank-one update of the cumulative statistics, both methods target the exact posterior distribution, from which a parameter draw can be obtained with every single observation. Simulation studies demonstrate the statistical and computational efficiency of the marginalized sampler. The methods are also applied to a real-world massive dataset on the incidentally truncated mortgage rates.

A Pólya–Gamma sampler for a generalized logistic regression

In this paper, we introduce a novel Bayesian data augmentation approach for estimating the parameters of the generalized logistic regression model. We propose a Pólya–Gamma sampler algorithm that allows us to sample from the exact posterior distribution, rather than relying on approximations. A simulation study illustrates the flexibility and accuracy of the proposed approach to capture heavy and light tails in binary response data of different dimensions. The algorithm performance is tested on simulated data. Furthermore, the methodology is applied to two different real datasets, where we demonstrate that the Pólya–Gamma sampler provides more precise estimates than the empirical likelihood method, outperforming approximate approaches.

Data-free likelihood-informed dimension reduction of Bayesian inverse problems

Identifying a low-dimensional informed parameter subspace offers a viable path to alleviating the dimensionality challenge in the sampled-based solution to large-scale Bayesian inverse problems. This paper introduces a novel gradient-based dimension reduction method in which the informed subspace does not depend on the data. This permits online–offline computational strategy where the expensive low-dimensional structure of the problem is detected in an offline phase, meaning before observing the data. This strategy is particularly relevant for multiple inversion problems as the same informed subspace can be reused. The proposed approach allows to control the approximation error (in expectation over the data) of the posterior distribution. We also present sampling strategies which exploit the informed subspace to draw efficiently samples from the exact posterior distribution. The method is successfully illustrated on two numerical examples: a PDE-based inverse problem with a Gaussian process prior and a tomography problem with Poisson data and a Besov- prior.

Exact posterior computation for the binomial–Kumaraswamy model

In Bayesian analysis, the well-known beta–binomial model is largely used as a conjugate structure, and the beta prior distribution is a natural choice to model parameters defined in the (0,1) range. The Kumaraswamy distribution has been used as a natural alternative to the beta distribution and has received great attention in statistics in the past few years, mainly due to the simplicity and the great variety of forms it can assume. However, the binomial–Kumaraswamy model is not conjugate, which may limit its use in situations where conjugacy is desired. This work provides the exact posterior distribution for the binomial–Kumaraswamy model using special functions. Besides the exact forms of the posterior moments, the predictive and the cumulative posterior distributions are provided. An example is used to illustrate the theory, in which the exact computation and the MCMC method are compared.

Using information theory to optimise epidemic models for real-time prediction and estimation.

The effective reproduction number, Rt, is a key time-varying prognostic for the growth rate of any infectious disease epidemic. Significant changes in Rt can forewarn about new transmissions within a population or predict the efficacy of interventions. Inferring Rt reliably and in real-time from observed time-series of infected (demographic) data is an important problem in population dynamics. The renewal or branching process model is a popular solution that has been applied to Ebola and Zika virus disease outbreaks, among others, and is currently being used to investigate the ongoing COVID-19 pandemic. This model estimates Rt using a heuristically chosen piecewise function. While this facilitates real-time detection of statistically significant Rt changes, inference is highly sensitive to the function choice. Improperly chosen piecewise models might ignore meaningful changes or over-interpret noise-induced ones, yet produce visually reasonable estimates. No principled piecewise selection scheme exists. We develop a practical yet rigorous scheme using the accumulated prediction error (APE) metric from information theory, which deems the model capable of describing the observed data using the fewest bits as most justified. We derive exact posterior prediction distributions for infected population size and integrate these within an APE framework to obtain an exact and reliable method for identifying the piecewise function best supported by available epidemic data. We find that this choice optimises short-term prediction accuracy and can rapidly detect salient fluctuations in Rt, and hence the infected population growth rate, in real-time over the course of an unfolding epidemic. Moreover, we emphasise the need for formal selection by exposing how common heuristic choices, which seem sensible, can be misleading. Our APE-based method is easily computed and broadly applicable to statistically similar models found in phylogenetics and macroevolution, for example. Our results explore the relationships among estimate precision, forecast reliability and model complexity.

Making Recursive Bayesian Inference Accessible

Bayesian models provide recursive inference naturally because they can formally reconcile new data and existing scientific information. However, popular use of Bayesian methods often avoids priors that are based on exact posterior distributions resulting from former studies. Two existing Recursive Bayesian methods are: Prior- and Proposal-Recursive Bayes. Prior-Recursive Bayes uses Bayesian updating, fitting models to partitions of data sequentially, and provides a way to accommodate new data as they become available using the posterior from the previous stage as the prior in the new stage based on the latest data. Proposal-Recursive Bayes is intended for use with hierarchical Bayesian models and uses a set of transient priors in first stage independent analyses of the data partitions. The second stage of Proposal-Recursive Bayes uses the posteriors from the first stage as proposals in a Markov chain Monte Carlo algorithm to fit the full model. We combine Prior- and Proposal-Recursive concepts to fit any Bayesian model, and often with computational improvements. We demonstrate our method with two case studies. Our approach has implications for big data, streaming data, and optimal adaptive design situations.

Bayesian Inference for ARFIMA Models

This article develops practical methods for Bayesian inference in the autoregressive fractionally integrated moving average (ARFIMA) model using the exact likelihood function, any proper prior distribution, and time series that may have thousands of observations. These methods utilize sequentially adaptive Bayesian learning, a sequential Monte Carlo algorithm that can exploit massively parallel desktop computing with graphics processing units (GPUs). The article identifies and solves several problems in the computation of the likelihood function that apparently have not been addressed in the literature. Four applications illustrate the utility of the approach. The most ambitious is an ARFIMA(2,d,2) model for the Campito tree ring time series (length 5405), for which the methods developed in the article provide an essentially uncorrelated sample of size 16,384 from the exact posterior distribution in under four hours. Less ambitious applications take as little as 4 minutes without exploiting GPUs.

Extending the Behrens-Fisher Problem to Testing Equality of Slopes in Linear Regression: The Bayesian Approach

Testing the equality of means of two normally distributed random variables when their variances are unequal is known in the statistical literature as the “Behrens-Fisher problem”. It is well-known that the posterior distributions of the parameters of interest are the primitive of Bayesian statistical inference. For routine implementation of statistical procedures based on posterior distributions, simple and efficient approaches are required. Since the computation of the exact posterior distribution of the Behrens-Fisher problem is obtained using numerical integration, several approximations are discussed and compared. Tests and Bayesian Highest-Posterior Density (H.P.D) intervals based upon these approximations are discussed. We extend the proposed approximations to test of parallelism in simple linear regression models.

A pruned recursive solution to the multiple change point problem

Long time series are often heterogeneous in nature. As such, the most appropriate model is one whose parameters are allowed to change through time. The exponential number of solutions to the multiple change point problem requires an efficient algorithm in order to be computationally feasible. Exact Bayesian solutions have at best quadratic complexity in the number of observations, which still can be too slow for very large data sets. Here, a pruned dynamic programming algorithm is proposed to fit a piecewise regression model with unknown break points to a data set. The algorithm removes unessential calculations, reducing the complexity of the most time consuming step of the algorithm from quadratic in the number of observations to quadratic in the average distance between change points. A distance measure is introduced that can be used to determine the divergence of the approximate joint posterior distribution from the exact posterior distribution. Analysis of two real data sets shows that this approximate algorithm produces a nearly identical representation of the joint posterior distribution on the locations of the change points, but with a significantly faster run time than its exact counterpart.