Marginal Likelihood Estimation Research Articles

Stabilizing training of affine coupling layers for high-dimensional variational inference

Abstract Variational inference with normalizing flows is an increasingly popular alternative to MCMC methods. In particular, normalizing flows based on affine coupling layers (Real NVPs) are frequently used due to their good empirical performance. In theory, increasing the depth of normalizing flows should lead to more accurate posterior approximations. However, in practice, training deep normalizing flows for approximating high-dimensional posterior distributions is often infeasible due to the high variance of the stochastic gradients. 
In this work, we show that previous methods for stabilizing the variance of stochastic gradient descent can be insufficient to achieve stable training of Real NVPs. As the source of the problem, we identify that, during training, samples often exhibit unusual high values. 
As a remedy, we propose a combination of two methods: (1) soft-thresholding of the scale in Real NVPs, and (2) a bijective soft log transformation of the samples. We evaluate these and other previously proposed modification on several challenging target distributions, including a high-dimensional horseshoe logistic regression model. Our experiments show that with our modifications, stable training of Real NVPs for posteriors with several thousand dimensions and heavy tails is possible, allowing for more accurate marginal likelihood estimation via importance sampling. 
Moreover, we evaluate several common training techniques and architecture choices and provide practical advise for training Real NVPs for high-dimensional variational inference. Finally, we also provide new empirical and theoretical justification that optimizing the evidence lower bound of normalizing flows leads to good posterior distribution coverage.

On Bayesian estimation of a latent trait model defined by a rank-based likelihood

Maximum likelihood estimation (frequentist) and Bayesian estimation are two common parameter estimation methods. However, maximum likelihood estimation faces limitations, including the effect of outliers, computational complexity, and issues with ordinal categorical data, leading to biased estimates and inaccurate coverage probabilities. To address these limitations, this study employed a latent trait model with a Bayesian marginal likelihood of rank-based estimation for parameter estimation. The simulation results demonstrated favorable performance of the proposed method. Trace plots of all parameters showed good distribution and quick convergence, with the potential scale reduction factor not exceeding 1, indicating no convergence issues. Furthermore, the posterior predictive check showed the simulated data closely resembled the observed data, indicating the method effectively captures within-region variation through a latent trait parameter. Performance metrics like mean absolute error, root mean square error, and 95% confidence interval coverage probability revealed the estimates from the proposed Bayesian method surpassed those from classical approaches. In conclusion, a latent traits model with Bayesian marginal likelihood and rank-based estimation is considered a superior parameter estimation technique compared to classical methods, particularly for dealing with ordinal categorical data.

Ab initio uncertainty quantification in scattering analysis of microscopy.

Estimating parameters from data is a fundamental problem in physics, customarily done by minimizing a loss function between a model and observed statistics. In scattering-based analysis, it is common to work in the reciprocal space. Researchers often employ their domain expertise to select a specific range of wave vectors for analysis, a choice that can vary depending on the specific case. We introduce another paradigm that defines a probabilistic generative model from the beginning of data processing and propagates the uncertainty for parameter estimation, termed the ab initio uncertainty quantification (AIUQ). As an illustrative example, we demonstrate this approach with differential dynamic microscopy (DDM) that extracts dynamical information through minimizing a loss function for the squared differences of the Fourier-transformed intensities, at a selected range of wave vectors. We first show that the conventional way of estimation in DDM is equivalent to fitting a temporal variogram in the reciprocal space using a latent factor model as the generative model. Then we derive the maximum marginal likelihood estimator, which optimally weighs the information at all wave vectors, therefore eliminating the need to select the range of wave vectors. Furthermore, we substantially reduce the computational cost of computing the likelihood function without approximation, by utilizing the generalized Schur algorithm for Toeplitz covariances. Simulated studies of a wide range of dynamical systems validate that the AIUQ method improves estimation accuracy and enables model selection with automated analysis. The utility of AIUQ is also demonstrated by three distinct sets of experiments: first in an isotropic Newtonian fluid, pushing limits of optically dense systems compared to multiple particle tracking; next in a system undergoing a sol-gel transition, automating the determination of gelling points and critical exponent; and lastly, in discerning anisotropic diffusive behavior of colloids in a liquid crystal. These studies demonstrate that the new approach does not require manually selecting the wave vector range and enables automated analysis.

Sequential and adaptive probabilistic integration for Bayesian model updating

The present study attempts to introduce a Bayesian perspective to marginal likelihood estimation. For this purpose, the Bayesian marginal likelihood inference (BMLI) framework is first developed, where the posterior statistics of the marginal likelihood are estimated using efficient random process sampling. The sequential and adaptive probabilistic integration (SAPI) method is then proposed within the BMLI framework to efficiently infer the posterior distribution of the model parameters as a result of the marginal likelihood estimation. In this method, an adaptive tempering of the likelihood function is presented to draw samples from a series of intermediate distributions to perform probabilistic integration sequentially, and the samples are finally converged to the target posterior domain. An active learning process is also designed to accelerate the sequential integration with a minimum number of likelihood function observations. Four numerical examples are investigated and the superiority of the SAPI method over other sequential sampling methods is demonstrated.

Bayesian model selection for structural damage identification: comparative analysis of marginal likelihood estimators

Bayesian model selection for structural damage identification: comparative analysis of marginal likelihood estimators

The Rank-2PL IRT Models for Forced-Choice Questionnaires: Maximum Marginal Likelihood Estimation with an EM Algorithm

The rank two-parameter logistic (Rank-2PL) item response theory models refer to a set of models applying the 2PL model in a sequential ranking process that occurs in forced-choice questionnaires. The multi-unidimensional pairwise preference with 2PL model (MUPP-2PL) is a Rank-2PL model for items with two statements. Focusing on items with three statements, we develop a maximum marginal likelihood estimation with an expectation-maximization algorithm to estimate item parameters and their standard errors. A simulation study is conducted to check parameter recovery, and then the model is applied to a real dataset. Finally, the findings are summarized and discussed, and future research is suggested.

Estimation Method of Covariance Matrix in Atmospheric Inversion of CO2 Emissions

Atmospheric inversion of CO2 Emissions is based on the correction of prior carbon dioxide flux inventories using concentration monitoring data and atmospheric transport models to obtain posterior carbon dioxide flux. In atmospheric inversion studies, fixed covariance functions are commonly used to generate covariance matrices, and the hyperparameters in the covariance functions are empirically estimated. In this study, we design and implement an ideal experiment based on meteorological data from the central urban area of Zhengzhou, using WRF-STILT to generate sensitivity matrices and construct real carbon emission inventories and prior inventories. Based on the real carbon emission inventories and sensitivity matrices of monitoring stations, simulated observation concentration values are generated. Firstly, based on the observed concentration values, sensitivity matrices of monitoring stations, prior inventories, and constructed covariance matrices, the values of hyperparameters are determined based on maximum marginal likelihood estimation. Then, the influence of different prior covariance functions on the inversion results is tested, and it is found that the prior covariance matrix generated by the balgovind covariance function is most suitable for the experimental data.

Detecting Episodic Evolution through Bayesian Inference of Molecular Clock Models.

Molecular evolutionary rate variation is a key aspect of the evolution of many organisms that can be modeled using molecular clock models. For example, fixed local clocks revealed the role of episodic evolution in the emergence of SARS-CoV-2 variants of concern. Like all statistical models, however, the reliability of such inferences is contingent on an assessment of statistical evidence. We present a novel Bayesian phylogenetic approach for detecting episodic evolution. It consists of computing Bayes factors, as the ratio of posterior and prior odds of evolutionary rate increases, effectively quantifying support for the effect size. We conducted an extensive simulation study to illustrate the power of this method and benchmarked it to formal model comparison of a range of molecular clock models using (log) marginal likelihood estimation, and to inference under a random local clock model. Quantifying support for the effect size has higher sensitivity than formal model testing and is straight-forward to compute, because it only needs samples from the posterior and prior distribution. However, formal model testing has the advantage of accommodating a wide range molecular clock models. We also assessed the ability of an automated approach, known as the random local clock, where branches under episodic evolution may be detected without their a priori definition. In an empirical analysis of a data set of SARS-CoV-2 genomes, we find "very strong" evidence for episodic evolution. Our results provide guidelines and practical methods for Bayesian detection of episodic evolution, as well as avenues for further research into this phenomenon.

Bayesian inverse problems with heterogeneous variance

AbstractWe consider inverse problems in Hilbert spaces under correlated Gaussian noise, and use a Bayesian approach to find their regularized solution. We focus on mildly ill‐posed inverse problems with fractional noise, using a novel wavelet‐based vaguelette–vaguelette approach. It allows us to apply sequence space methods without assuming that all operators are simultaneously diagonalizable. The results are proved for more general bases and covariance operators. Our primary aim is to study posterior contraction rate in such inverse problems over Sobolev classes and compare it to the derived minimax rate. Secondly, we study effect of plugging in a consistent estimator of variances in sequence space on the posterior contraction rate. This result is applied to the problem with error in forward operator. Thirdly, we show that empirical Bayes posterior distribution with a plugged‐in maximum marginal likelihood estimator of the prior scale contracts at the optimal rate, adaptively, in the minimax sense.

A double-Cox model for non-proportional hazards survival analysis with frailty.

The Cox regression, a semi-parametric method of survival analysis, is extremely popular in biomedical applications. The proportional hazards assumption is a key requirement in the Cox model. To accommodate non-proportional hazards, we propose to parameterize the shape parameter of the baseline hazard function using the additional, separate Cox-regression term which depends on the vector of the covariates. This parametrization retains the general form of the hazard function over the strata and is similar to one in Devarajan and Ebrahimi (Comput Stat Data Anal. 2011;55:667-676) in the case of the Weibull distribution, but differs for other hazard functions. We call this model the double-Cox model. We formally introduce the double-Cox model with shared frailty and investigate, by simulation, the estimation bias and the coverage of the proposed point and interval estimation methods for the Gompertz and the Weibull baseline hazards. For real-life applications with low frailty variance and a large number of clusters, the marginal likelihood estimation is almost unbiased and the profile likelihood-based confidence intervals provide good coverage for all model parameters. We also compare the results from the over-parametrized double-Cox model to those from the standard Cox model with frailty in the case of the scale-only proportional hazards. The model is illustrated on an example of the survival after a diagnosis of type 2 diabetes mellitus. The R programs for fitting the double-Cox model are available on Github.

Target-aware Bayesian inference via generalized thermodynamic integration

In Bayesian inference, we are usually interested in the numerical approximation of integrals that are posterior expectations or marginal likelihoods (a.k.a., Bayesian evidence). In this paper, we focus on the computation of the posterior expectation of a function f(textbf{x}). We consider a target-aware scenario where f(textbf{x}) is known in advance and can be exploited in order to improve the estimation of the posterior expectation. In this scenario, this task can be reduced to perform several independent marginal likelihood estimation tasks. The idea of using a path of tempered posterior distributions has been widely applied in the literature for the computation of marginal likelihoods. Thermodynamic integration, path sampling and annealing importance sampling are well-known examples of algorithms belonging to this family of methods. In this work, we introduce a generalized thermodynamic integration (GTI) scheme which is able to perform a target-aware Bayesian inference, i.e., GTI can approximate the posterior expectation of a given function. Several scenarios of application of GTI are discussed and different numerical simulations are provided.

Comparing the transmission potential from sequence and surveillance data of 2009 North American influenza pandemic waves.

Technological advancements in phylodynamic modeling coupled with the accessibility of real-time pathogen genetic data are increasingly important for understanding the infectious disease transmission dynamics. In this study, we compare the transmission potentials of North American influenza A(H1N1)pdm09 derived from sequence data to that derived from surveillance data. The impact of the choice of tree-priors, informative epidemiological priors, and evolutionary parameters on the transmission potential estimation is evaluated. North American Influenza A(H1N1)pdm09 hemagglutinin (HA) gene sequences are analyzed using the coalescent and birth-death tree prior models to estimate the basic reproduction number (R 0 ). Epidemiological priors gathered from published literature are used to simulate the birth-death skyline models. Path-sampling marginal likelihood estimation is conducted to assess model fit. A bibliographic search to gather surveillance-based R 0 values were consistently lower (mean≤1.2) when estimated by coalescent models than by the birth-death models with informative priors on the duration of infectiousness (mean≥1.3 to ≤2.88 days). The user-defined informative priors for use in the birth-death model shift the directionality of epidemiological and evolutionary parameters compared to non-informative estimates. While there was no certain impact of clock rate and tree height on the R 0 estimation, an opposite relationship was observed between coalescent and birth-death tree priors. There was no significant difference (p=0.46) between the birth-death model and surveillance R 0 estimates. This study concludes that tree-prior methodological differences may have a substantial impact on the transmission potential estimation as well as the evolutionary parameters. The study also reports a consensus between the sequence-based R 0 estimation and surveillance-based R 0 estimates. Altogether, these outcomes shed light on the potential role of phylodynamic modeling to augment existing surveillance and epidemiological activities to better assess and respond to emerging infectious diseases.

Bayesian optimization of hyperparameters from noisy marginal likelihood estimates

SummaryBayesian models often involve a small set of hyperparameters determined by maximizing the marginal likelihood. Bayesian optimization is an iterative method where a Gaussian process posterior of the underlying function is sequentially updated by new function evaluations. We propose a novel Bayesian optimization framework for situations where the user controls the computational effort and therefore the precision of the function evaluations. This is a common situation in econometrics where the marginal likelihood is often computed by Markov chain Monte Carlo or importance sampling methods. The new acquisition strategy gives the optimizer the option to explore the function with cheap noisy evaluations and therefore find the optimum faster. The method is applied to estimating the prior hyperparameters in two popular models on US macroeconomic time series data: the steady‐state Bayesian vector autoregressive (BVAR) and the time‐varying parameter BVAR with stochastic volatility.

Posterior contraction and testing for multivariate isotonic regression

We consider the nonparametric regression problem with multiple predictors and an additive error, where the regression function is assumed to be coordinatewise nondecreasing. We propose a Bayesian approach to make inferences on the multivariate monotone regression function, obtain the posterior contraction rate, and construct a universally consistent Bayesian testing procedure for multivariate monotonicity. To facilitate posterior analysis, we temporarily set aside the shape restrictions, and endow a prior on blockwise constant regression functions with independently normally distributed heights. The unknown variance of the error term is either estimated by the marginal maximum likelihood estimate, or is equipped with an inverse-gamma prior. Then the unrestricted block-heights are a posteriori also independently normally distributed given the error variance, by conjugacy. To comply with the shape restrictions, we project samples from the unrestricted posterior onto the class of multivariate monotone functions, inducing the “projection-posterior distribution”, to be used for making an inference. Under an L1-metric, we show that the projection-posterior based on n independent samples contracts around the true monotone regression function at the optimal rate n−1∕(2+d). Then we construct a Bayesian test for multivariate monotonicity based on the posterior probability of a shrinking neighborhood of the class of multivariate monotone functions. We show that the test is universally consistent, that is, the level of the Bayesian test goes to zero, and the power at any fixed alternative goes to one. Moreover, we show that for a smooth alternative function, power goes to one as long as its distance to the class of multivariate monotone functions is at least of the order of the estimation error for a smooth function. To the best of our knowledge, no other test for multivariate monotonicity is available in the Bayesian or the frequentist literature.

Comparing stochastic volatility specifications for large Bayesian VARs

Large Bayesian vector autoregressions with various forms of stochastic volatility have become increasingly popular in empirical macroeconomics. One main difficulty for practitioners is to choose the most suitable stochastic volatility specification for their particular application. We develop Bayesian model comparison methods–based on marginal likelihood estimators that combine conditional Monte Carlo and adaptive importance sampling–to choose among a variety of stochastic volatility specifications. The proposed methods can also be used to select an appropriate shrinkage prior on the VAR coefficients, which is a critical component for avoiding over-fitting in high-dimensional settings. Using US quarterly data of different dimensions, we find that both the Cholesky stochastic volatility and factor stochastic volatility outperform the common stochastic volatility specification. Their superior performance, however, can mostly be attributed to the more flexible priors that accommodate cross-variable shrinkage.

Bayesian updating and marginal likelihood estimation by cross entropy based importance sampling

We introduce a novel simulation-based method for Bayesian analysis to learn model parameters based on data. The method employs importance sampling (IS) to construct a sample-based approximation of the posterior probability density function (PDF) and estimate the marginal likelihood. We propose to build the IS density through an adaptive sampling approach based on the cross entropy (CE) method. The aim is to identify the parameters of a chosen parametric distribution family that minimize its Kullback-Leibler divergence from the target posterior PDF. An adaptive multi-level approach, based on tempering of the likelihood function, is proposed to efficiently solve this CE optimization problem. We investigate the appropriate choice of the parametric distribution, depending on the number of uncertain model parameters and nature of the posterior density. Numerical studies demonstrate the performance of the proposed method.

A general class of small area estimation using calibrated hierarchical likelihood approach with applications to COVID-19 data

The direct estimation techniques in small area estimation (SAE) models require sufficiently large sample sizes to provide accurate estimates. Hence, indirect model-based methodologies are developed to incorporate auxiliary information. The most commonly used SAE models, including the Fay-Herriot (FH) model and its extended models, are estimated using marginal likelihood estimation and the Bayesian methods, which rely heavily on the computationally intensive integration of likelihood function. In this article, we propose a Calibrated Hierarchical (CH) likelihood approach to obtain SAE through hierarchical estimation of fixed effects and random effects with the regression calibration method for bias correction. The latent random variables at the domain level are treated as ‘parameters’ and estimated jointly with other parameters of interest. Then the dispersion parameters are estimated iteratively based on the Laplace approximation of the profile likelihood. The proposed method avoids the intractable integration to estimate the marginal distribution. Hence, it can be applied to a wide class of distributions, including generalized linear mixed models, survival analysis, and joint modeling with distinct distributions. We demonstrate our method using an area-level analysis of publicly available count data from the novel coronavirus (COVID-19) positive cases.

Development of Integrated Choice and Latent Variable (ICLV) Models Using Matrix-Based Analytic Approximation and Automatic Differentiation Methods on TensorFlow Platform

This study further explores the multinomial probit-based integrated choice and latent variable (ICLV) models. The LDLT matrix-based analytic approximation methods, including Mendell and Elston (ME) method, bivariate ME (BME) method, and two-variate bivariate screening (TVBS) method, were adapted to calculate the multivariate cumulative normal distribution (MVNCD) function in the ICLV model because of the better performances in accuracy and computational time. Integrated with the composite marginal likelihood (CML) estimation approach, the ICLV model based on high-dimensional integration can be estimated accurately within a reasonable time. In this study, some three-alternative and four-alternative ICLV models are simulated to examine their abilities to recover model parameters. It is found that the parameter estimates and standard error estimates are acceptable for both models and the computational time is expected to decrease using tensor data structures on the TensorFlow platform. For the four-alternative ICLV models, the TVBS method has the highest level of accuracy. The BME method is also a good alternative to TVBS if computational time is of great concern. The application of the automatic differentiation (AD) technique in the model can free researchers from coding analytical gradients of log-likelihood functions and thereby greatly reduce the workload of researchers.

The node-wise Pseudo-marginal method: model selection with spatial dependence on latent graphs

Motivated by problems from neuroimaging in which existing approaches make use of “mass univariate” analysis which neglects spatial structure entirely, but the full joint modelling of all quantities of interest is computationally infeasible, a novel method for incorporating spatial dependence within a (potentially large) family of model-selection problems is presented. Spatial dependence is encoded via a Markov random field model for which a variant of the pseudo-marginal Markov chain Monte Carlo algorithm is developed and extended by a further augmentation of the underlying state space. This approach allows the exploitation of existing unbiased marginal likelihood estimators used in settings in which spatial independence is normally assumed thereby facilitating the incorporation of spatial dependence using non-spatial estimates with minimal additional development effort. The proposed algorithm can be realistically used for analysis of moderately sized data sets such as 2D slices of whole 3D dynamic PET brain images or other regions of interest. Principled approximations of the proposed method, together with simple extensions based on the augmented spaces, are investigated and shown to provide similar results to the full pseudo-marginal method. Such approximations and extensions allow the improved performance obtained by incorporating spatial dependence to be obtained at negligible additional cost. An application to measured PET image data shows notable improvements in revealing underlying spatial structure when compared to current methods that assume spatial independence.

Inference and model determination for temperature-driven non-linear ecological models

This paper is concerned with a contemporary Bayesian approach to the effect of temperature on developmental rates. We develop statistical methods using recent computational tools to model four commonly used ecological non-linear mathematical curves that describe arthropods’ developmental rates. Such models address the effect of temperature fluctuations on the developmental rate of arthropods. In addition to the widely used Gaussian distributional assumption, we also explore Inverse Gamma-based alternatives, which naturally accommodate adaptive variance fluctuation with temperature. Moreover, to overcome the associated parameter indeterminacy in the case of no development, we suggest the zero-inflated Inverse Gamma model. The ecological models are compared graphically via posterior predictive plots and quantitatively via marginal likelihood estimates and Information criteria. Inference is performed using the Stan software and we investigate the statistical and computational efficiency of its Hamiltonian Monte Carlo and Variational Inference methods. We also explore model uncertainty and employ Bayesian Model Averaging framework for robust estimation of the key ecological parameters.