An adaptive Kriging surrogate method for efficient joint estimation of hydraulic and biochemical parameters in reactive transport modeling

The validity of using Gaussian assumptions for model residuals in uncertainty quantification of a groundwater reactive transport model was evaluated in this study. Least squares regression methods explicitly assume Gaussian residuals, and the assumption leads to Gaussian likelihood functions, model parameters, and model predictions. While the Bayesian methods do not explicitly require the Gaussian assumption, Gaussian residuals are widely used. This paper shows that the residuals of the reactive transport model are non-Gaussian, heteroscedastic, and correlated in time; characterizing them requires using a generalized likelihood function such as the formal generalized likelihood function developed by Schoups and Vrugt (2010). For the surface complexation model considered in this study for simulating uranium reactive transport in groundwater, parametric uncertainty is quantified using the least squares regression methods and Bayesian methods with both Gaussian and formal generalized likelihood functions. While the least squares methods and Bayesian methods with Gaussian likelihood function produce similar Gaussian parameter distributions, the parameter distributions of Bayesian uncertainty quantification using the formal generalized likelihood function are non-Gaussian. In addition, predictive performance of formal generalized likelihood function is superior to that of least squares regression and Bayesian methods with Gaussian likelihood function. The Bayesian uncertainty quantification is conducted using the differential evolution adaptive metropolis (DREAM(zs)) algorithm; as a Markov chain Monte Carlo (MCMC) method, it is a robust tool for quantifying uncertainty in groundwater reactive transport models. For the surface complexation model, the regression-based local sensitivity analysis and Morris- and DREAM(ZS)-based global sensitivity analysis yield almost identical ranking of parameter importance. The uncertainty analysis may help select appropriate likelihood functions, improve model calibration, and reduce predictive uncertainty in other groundwater reactive transport and environmental modeling.

In diffusion MRI analysis, advances in biophysical multi-compartment modeling have gained popularity over the conventional Diffusion Tensor Imaging (DTI), because they can obtain a greater specificity in relating the dMRI signal to underlying cellular microstructure. Biophysical multi-compartment models require a parameter estimation, typically performed using either the Maximum Likelihood Estimation (MLE) or the Markov Chain Monte Carlo (MCMC) sampling. Whereas, the MLE provides only a point estimate of the fitted model parameters, the MCMC recovers the entire posterior distribution of the model parameters given in the data, providing additional information such as parameter uncertainty and correlations. MCMC sampling is currently not routinely applied in dMRI microstructure modeling, as it requires adjustment and tuning, specific to each model, particularly in the choice of proposal distributions, burn-in length, thinning, and the number of samples to store. In addition, sampling often takes at least an order of magnitude, more time than non-linear optimization. Here we investigate the performance of the MCMC algorithm variations over multiple popular diffusion microstructure models, to examine whether a single, well performing variation could be applied efficiently and robustly to many models. Using an efficient GPU-based implementation, we showed that run times can be removed as a prohibitive constraint for the sampling of diffusion multi-compartment models. Using this implementation, we investigated the effectiveness of different adaptive MCMC algorithms, burn-in, initialization, and thinning. Finally we applied the theory of the Effective Sample Size, to the diffusion multi-compartment models, as a way of determining a relatively general target for the number of samples needed to characterize parameter distributions for different models and data sets. We conclude that adaptive Metropolis methods increase MCMC performance and select the Adaptive Metropolis-Within-Gibbs (AMWG) algorithm as the primary method. We furthermore advise to initialize the sampling with an MLE point estimate, in which case 100 to 200 samples are sufficient as a burn-in. Finally, we advise against thinning in most use-cases and as a relatively general target for the number of samples, we recommend a multivariate Effective Sample Size of 2,200.

Markov Chain Monte Carlo (MCMC) has been essential in tracking vehicle undergoing disturbances for traffic surveillance purposes. It is capable of tracking vehicle by estimating the vehicle's position with the sampling of probability distributions. However the accuracy of the position estimation is highly dependent on the sampling efficiency of MCMC. Therefore the sample size of the MCMC is adapted to track the vehicle according to the disturbances encountered. The adaptive sample size of MCMC is determined by using the CUSUM path plot and variance ratio convergence diagnostic algorithm. To further enhance the convergence speed, genetic crossover and mutation operator is introduced into the adaptive MCMC. The genetic operator (GO) is capable of reduces the variance between samples and hence allowing faster convergence speed on the MCMC samples. Experimental results have shown that the GO adaptive MCMC tracking algorithm have better tracking performances with consumption of lesser sample size.

BackgroundOne of the goals of the Systems Biology community is to have a detailed map of all biological interactions in an organism. One small yet important step in this direction is the creation of biological networks from post-genomic data. Bayesian networks are a very promising model for the inference of regulatory networks in Systems Biology. Usually, Bayesian networks are sampled with a Markov Chain Monte Carlo (MCMC) sampler in the structure space. Unfortunately, conventional MCMC sampling schemes are often slow in mixing and convergence. To improve MCMC convergence, an alternative method is proposed and tested with different sets of data. Moreover, the proposed method is compared with the traditional MCMC sampling scheme.ResultsIn the proposed method, a simpler and faster method for the inference of regulatory networks, Graphical Gaussian Models (GGMs), is integrated into the Bayesian network inference, trough a Hierarchical Bayesian model. In this manner, information about the structure obtained from the data with GGMs is taken into account in the MCMC scheme, thus improving mixing and convergence. The proposed method is tested with three types of data, two from simulated models and one from real data. The results are compared with the results of the traditional MCMC sampling scheme in terms of network recovery accuracy and convergence. The results show that when compared with a traditional MCMC scheme, the proposed method presents improved convergence leading to better network reconstruction with less MCMC iterations.ConclusionsThe proposed method is a viable alternative to improve mixing and convergence of traditional MCMC schemes. It allows the use of Bayesian networks with an MCMC sampler with less iterations. The proposed method has always converged earlier than the traditional MCMC scheme. We observe an improvement in accuracy of the recovered networks for the Gaussian simulated data, but this improvement is absent for both real data and data simulated from ODE.Electronic supplementary materialThe online version of this article (doi:10.1186/s12859-015-0734-6) contains supplementary material, which is available to authorized users.

Environmental protection is one of the themes of today's world. The forest is a recycler of carbon dioxide and natural oxygen bar. Protection of forests, monitoring of forest growth is long-term task of environmental protection. It is very important to automatically statistic the forest coverage rate using optical remote sensing images and the computer, by which we can timely understand the status of the forest of an area, and can be freed from tedious manual statistics. Towards the problem of computational complexity of the global optimization using convexification, this paper proposes a level set segmentation method based on Markov chain Monte Carlo (MCMC) sampling and applies it to forest segmentation in remote sensing images. The presented method needs not to do any convexity transformation for the energy functional of the goal, and uses MCMC sampling method with global optimization capability instead. The possible local minima occurring by using gradient descent method is also avoided. There are three major contributions in the paper. Firstly, by using MCMC sampling, the convexity of the energy functional is no longer necessary and global optimization can still be achieved. Secondly, taking advantage of the data (texture) and knowledge (a priori color) to guide the construction of Markov chain, the convergence rate of Markov chains is improved significantly. Finally, the level set segmentation method by integrating a priori color and texture for forest is proposed. The experiments show that our method can efficiently and accurately segment forest in remote sensing images.

Bayesian inference under a set of priors, called robust Bayesian analysis, allows for estimation of parameters within a model and quantification of epistemic uncertainty in quantities of interest by bounded (or imprecise) probability. Iterative importance sampling can be used to estimate bounds on the quantity of interest by optimizing over the set of priors. A method for iterative importance sampling when the robust Bayesian inference relies on Markov chain Monte Carlo (MCMC) sampling is proposed. To accommodate the MCMC sampling in iterative importance sampling, a new expression for the effective sample size of the importance sampling is derived, which accounts for the correlation in the MCMC samples. To illustrate the proposed method for robust Bayesian analysis, iterative importance sampling with MCMC sampling is applied to estimate the lower bound of the overall effect in a previously published meta-analysis with a random effects model. The performance of the method compared to a grid search method and under different degrees of prior-data conflict is also explored.

Sparsity has become a key concept for solving of high-dimensional inverse problems using variational regularization techniques. Recently, using similar sparsity-constraints in the Bayesian framework for inverse problems by encoding them in the prior distribution has attracted attention. Important questions about the relation between regularization theory and Bayesian inference still need to be addressed when using sparsity promoting inversion. A practical obstacle for these examinations is the lack of fast posterior sampling algorithms for sparse, high-dimensional Bayesian inversion. Accessing the full range of Bayesian inference methods requires being able to draw samples from the posterior probability distribution in a fast and efficient way. This is usually done using Markov chain Monte Carlo (MCMC) sampling algorithms. In this paper, we develop and examine a new implementation of a single component Gibbs MCMC sampler for sparse priors relying on L1-norms. We demonstrate that the efficiency of our Gibbs sampler increases when the level of sparsity or the dimension of the unknowns is increased. This property is contrary to the properties of the most commonly applied Metropolis–Hastings (MH) sampling schemes. We demonstrate that the efficiency of MH schemes for L1-type priors dramatically decreases when the level of sparsity or the dimension of the unknowns is increased. Practically, Bayesian inversion for L1-type priors using MH samplers is not feasible at all. As this is commonly believed to be an intrinsic feature of MCMC sampling, the performance of our Gibbs sampler also challenges common beliefs about the applicability of sample based Bayesian inference.

We address the problem of parameter estimation of chirplets which are chirp signals with Gaussian shaped envelopes. The procedure we propose is an extension of our previous work on estimation of chirp signals (Lin and Djuric, 2000), and it is based on MCMC sampling. For fast convergence of the Markov chain Monte Carlo (MCMC) sampling based method, a critical step is the initialization of the method Since the chirplets have finite durations and may or may not overlap in time, we propose initialization procedures for each of these cases. We have tested the method by extensive simulations and compared it with Cramer-Rao bounds. The obtained results have been excellent.

Summary Estimating abundance from mark–recapture data is challenging when capture probabilities vary among individuals. Initial solutions to this problem were based on fitting conditional likelihoods and estimating abundance as a derived parameter. More recently, Bayesian methods using full likelihoods have been implemented via reversible jump Markov chain Monte Carlo sampling (RJMCMC) or data augmentation (DA). The latter approach is easily implemented in available software and has been applied to fit models that allow for heterogeneity in both open and closed populations. However, both RJMCMC and DA may be inefficient when modelling large populations. We describe an alternative approach using Monte Carlo (MC) integration to approximate the posterior density within a Markov chain Monte Carlo (MCMC) sampling scheme. We show how this Monte Carlo within MCMC (MCWM) approach may be used to fit a simple, closed population model including a single individual covariate and present results from a simulation study comparing RJMCMC, DA and MCWM. We found that MCWM can provide accurate inference about population size and can be more efficient than both RJMCMC and DA. The efficiency of MCWM can also be improved by using advanced MC methods like antithetic sampling. Finally, we apply MCWM to estimate the abundance of meadow voles (Microtus pennsylvanicus) at the Patuxent Wildlife Research Center in 1982 allowing for capture probabilities to vary as a function body mass.

Complex-valued Bayesian parameter estimation via Markov chain Monte Carlo

[1] Bayesian analysis has become vital to uncertainty quantification in groundwater modeling, but its application has been hindered by the computational cost associated with numerous model executions required by exploring the posterior probability density function (PPDF) of model parameters. This is particularly the case when the PPDF is estimated using Markov Chain Monte Carlo (MCMC) sampling. In this study, a new approach is developed to improve the computational efficiency of Bayesian inference by constructing a surrogate of the PPDF, using an adaptive sparse-grid high-order stochastic collocation (aSG-hSC) method. Unlike previous works using first-order hierarchical basis, this paper utilizes a compactly supported higher-order hierarchical basis to construct the surrogate system, resulting in a significant reduction in the number of required model executions. In addition, using the hierarchical surplus as an error indicator allows locally adaptive refinement of sparse grids in the parameter space, which further improves computational efficiency. To efficiently build the surrogate system for the PPDF with multiple significant modes, optimization techniques are used to identify the modes, for which high-probability regions are defined and components of the aSG-hSC approximation are constructed. After the surrogate is determined, the PPDF can be evaluated by sampling the surrogate system directly without model execution, resulting in improved efficiency of the surrogate-based MCMC compared with conventional MCMC. The developed method is evaluated using two synthetic groundwater reactive transport models. The first example involves coupled linear reactions and demonstrates the accuracy of our high-order hierarchical basis approach in approximating high-dimensional posteriori distribution. The second example is highly nonlinear because of the reactions of uranium surface complexation, and demonstrates how the iterative aSG-hSC method is able to capture multimodal and non-Gaussian features of PPDF caused by model nonlinearity. Both experiments show that aSG-hSC is an effective and efficient tool for Bayesian inference.

This paper addresses the problem of estimating the lower atmospheric refractivity (M profile) under nonstandard propagation conditions frequently encountered in low‐altitude maritime radar applications. This is done by statistically estimating the duct strength (range‐ and height‐dependent atmospheric index of refraction) from the sea surface reflected radar clutter. These environmental statistics can then be used to predict the radar performance. In previous work, genetic algorithms (GA) and Markov chain Monte Carlo (MCMC) samplers were used to calculate the atmospheric refractivity from returned radar clutter. Although GA is fast and estimates the maximum a posteriori (MAP) solution well, it poorly calculates the multidimensional integrals required to obtain the means, variances, and underlying posterior probability distribution functions of the estimated parameters. More accurate distributions and integral calculations can be obtained using MCMC samplers, such as the Metropolis‐Hastings and Gibbs sampling (GS) algorithms. Their drawback is that they require a large number of samples relative to the global optimization techniques such as GA and become impractical with an increasing number of unknowns. A hybrid GA‐MCMC method based on the nearest neighborhood algorithm is implemented in this paper. It is an improved GA method which improves integral calculation accuracy through hybridization with a MCMC sampler. Since the number of forward models is determined by GA, it requires fewer forward model samples than a MCMC, enabling inversion of atmospheric models with a larger number of unknowns.

Global illumination algorithms using Markov chain Monte Carlo (MCMC) sampling are well-known for their efficiency in scenes with complex light transport. Samples in such algorithms are generated as a history of Markov chain states so that they are distributed according to the contributions to the image. The whole process is done based only on the information of the path contributions and user-defined transition probabilities from one state to the others. In light transport simulation, however, there is more information that can be used to improve the efficiency of path sampling. A notable example is multiple importance sampling (MIS) in bidirectional path tracing, which utilizes the probability densities of constructing a given path with different estimators. While MIS is a powerful ordinary Monte Carlo method, how to incorporate such additional information into MCMC sampling has been an open problem. We introduce a novel MCMC sampling framework, primary space serial tempering , which fuses the ideas of MCMC sampling and MIS for the first time. The key idea is to explore not only the sample space using a Markov chain, but also different estimators to generate samples by utilizing the information already available for MIS. Based on this framework, we also develop a novel rendering algorithm, multiplexed Metropolis light transport , which automatically and adaptively constructs paths with appropriate techniques as predicted by MIS. The final algorithm is very easy to implement, yet in many cases shows comparable (or even better) performance than significantly more complex MCMC rendering algorithms.

In Bayesian inference, a joint posterior distribution is available through the likelihood function and a prior distribution. One purpose of Bayesian inference is to calculate and display marginal posterior densities because the marginal posterior densities provide complete information about parameters of interest. As shown in Chapter 2, a Markov chain Monte Carlo (MCMC) sampling algorithm, such as the Gibbs sampler or a Metropolis-Hastings algorithm, can be used to draw MCMC samples from the posterior distribution. Chapter 3 also demonstrates how we can easily obtain posterior quantities such as posterior means, posterior standard deviations, and other posterior quantities from MCMC samples. However, when a Bayesian model becomes complicated, it may be difficult to obtain a reliable estimator of a marginal posterior density based on the MCMC sample. A traditional method for estimating marginal posterior densities is kernel density estimation. Since the kernel density estimator is nonparametric, it may not be efficient. On the other hand, the kernel density estimator may not be applicable for some complicated Bayesian models. In the context of Bayesian inference, the joint posterior density is typically known up to a normalizing constant. Using the structure of a posterior density, a number of authors (e.g., Gelfand, Smith, and Lee 1992; Johnson 1992; Chen 1993 and 1994; Chen and Shao 1997c; Chib 1995; Verdinelli and Wasserman 1995) propose parametric marginal posterior density estimators based on the MCMC sample. In this chapter, we present several available Monte Carlo (MC) methods for computing marginal posterior density estimators, and we also discuss how well marginal posterior density estimation works using the Kullback—Leibler (K—L) divergence as a performance measure.

I describe algorithms for drawing from distributions using adaptive Markov chain Monte Carlo (MCMC) methods; I introduce a Mata function for performing adaptive MCMC, amcmc(); and I present a suite of functions, amcmc_ *(), that allows an alternative implementation of adaptive MCMC. amcmc() and amcmc_ *() can be used with models set up to work with Mata's moptimize( ) (see [M-5] moptimize( )) or optimize( ) (see [M-5] optimize( )) or with standalone functions. To show how the routines can be used in estimation problems, I give two examples of what Chernozhukov and Hong (2003, Journal of Econometrics 115: 293–346) refer to as quasi-Bayesian or Laplace-type estimators—simulation-based estimators using MCMC sampling. In the first example, I illustrate basic ideas and show how a simple linear model can be fit by simulation. In the next example, I describe simulation-based estimation of a censored quantile regression model following Powell (1986, Journal of Econometrics 32: 143–155); the discussion describes the workings of the command mcmccqreg. I also present an example of how the routines can be used to draw from distributions without a normalizing constant and used in Bayesian estimation of a mixed logit model. This discussion introduces the command bayesmixedlogit.