Monte Carlo Sampling Algorithm Research Articles

Bayesian Joint Modeling of Single-Cell Expression Data and Bulk Spatial Transcriptomic Data

Single-cell RNA-sequencing (scRNA-seq) enables gene expression profiling at single-cell resolution, but it loses the spatial information of cells for solid tissues during the tissue dissociation step before sequencing. In contrast, bulk spatial transcriptomics (ST) methods can measure the expression of spatially organized spots in solid tissues, but as a spot comprises dozens of cells, ST expression levels are averaged signals and lack cellular resolution. Joint analysis of these two complementary data types provides the opportunity to recover the spatial patterns of cell types and obtain the cellular enrichment of spots. However, there is a lack of unified statistical methods to achieve this goal. This study develops a Bayesian statistical method named BEATS to jointly model scRNA-seq data and bulk ST data from a common sample in the presence of cellular and spatial heterogeneity. BEATS can simultaneously (a) discover cell types, where cells in a cell type share mean expression profiles; (b) identify spot regions, where a region is a set of spots with the same cellular compositions; and (c) estimate cell-type proportions for each spot region. The Bayesian posterior inference is performed through a hybrid Markov chain Monte Carlo sampling algorithm. Extensive simulation studies and application to datasets on pancreatic ductal adenocarcinoma tissues demonstrate the practical utility of BEATS.

Bayesian network meta-regression hierarchical models using heavy-tailed multivariate random effects with covariate-dependent variances.

Network meta-analysis (NMA) is gaining popularity in evidence synthesis and network meta-regression allows us to incorporate potentially important covariates into network meta-analysis. In this article, we propose a Bayesian network meta-regression hierarchical model and assume a general multivariate t distribution for the random treatment effects. The multivariate t distribution is desired for heavy-tailed random effects and converges to the multivariate normal distribution when the degrees of freedom go to infinity. Moreover, in NMA, some treatments are compared only in a single study. To overcome such sparsity, we propose a log-linear regression model for the variances of the random effects and incorporate aggregate covariates into modeling the variance components. We develop a Markov chain Monte Carlo sampling algorithm to sample from the posterior distribution via the collapsed Gibbs technique. We further use the deviance information criterion and the logarithm of the pseudo-marginal likelihood for model comparison. A simulation study is conducted and a detailed analysis from our motivating case study is carried out to further demonstrate the proposed methodology.

Bayesian Analysis of Intraday Stochastic Volatility Models of High-Frequency Stock Returns with Skew Heavy-Tailed Errors

Intraday high-frequency data of stock returns exhibit not only typical characteristics (e.g., volatility clustering and the leverage effect) but also a cyclical pattern of return volatility that is known as intraday seasonality. In this paper, we extend the stochastic volatility (SV) model for application with such intraday high-frequency data and develop an efficient Markov chain Monte Carlo (MCMC) sampling algorithm for Bayesian inference of the proposed model. Our modeling strategy is two-fold. First, we model the intraday seasonality of return volatility as a Bernstein polynomial and estimate it along with the stochastic volatility simultaneously. Second, we incorporate skewness and excess kurtosis of stock returns into the SV model by assuming that the error term follows a family of generalized hyperbolic distributions, including variance-gamma and Student’s t distributions. To improve efficiency of MCMC implementation, we apply an ancillarity-sufficiency interweaving strategy (ASIS) and generalized Gibbs sampling. As a demonstration of our new method, we estimate intraday SV models with 1 min return data of a stock price index (TOPIX) and conduct model selection among various specifications with the widely applicable information criterion (WAIC). The result shows that the SV model with the skew variance-gamma error is the best among the candidates.

Scalable Fully Bayesian Gaussian Process Modeling and Calibration With Adaptive Sequential Monte Carlo for Industrial Applications

Abstract Gaussian process (GP) regression or kriging has been extensively applied in the engineering literature for the purposes of building a cheap-to-evaluate surrogate, within the contexts of multi-fidelity modeling, model calibration, and design optimization. With the ongoing automation of manufacturing and industrial practices as a part of Industry 4.0, there has been a greater need for advancing GP regression techniques to handle challenges such as high input dimensionality, data paucity or big data problems, these consist primarily of proposing efficient design of experiments, optimal data acquisition strategies, sparsifying covariance kernels, and other mathematical tricks. In this work, our attention is focused on the challenges of efficiently training a GP model, which, to the authors opinion, has attracted very little attention and is to-date poorly addressed. The performance of widely used training approaches such as maximum likelihood estimation and Markov Chain Monte Carlo (MCMC) sampling can deteriorate significantly in high-dimensional and big data problems and can lead to cost deficient implementations of critical importance to many industrial applications. Here, we compare an Adaptive Sequential Monte Carlo (ASMC) sampling algorithm to classic MCMC sampling strategies and we demonstrate the effectiveness of our implementation on several mathematical problems and challenging industry applications of varying complexity. The computational time savings of the ASMC approach manifest in large-scale problems helping us to push the boundaries of applicability and scalability of GPs for model calibration in various domains of the industry, including but not limited to design automation, design engineering, smart manufacturing, predictive maintenance, and supply chain manufacturing.

Quantifying input error in hydrologic modeling using the Bayesian error analysis with reordering (BEAR) approach

Input error in hydrologic models, mainly arising from observed precipitation, can impair model calibration. Precise and reliable identification of input error is important for improved model parameter estimation. However, limited information about the nature of the error and the memory of the hydrologic system make it challenging to disentangle input error from the total residual error. Based on a Sequential Monte Carlo sampling algorithm, Bayesian error analysis with reordering (BEAR) method is developed to quantify the unknown input error series by sampling errors from an assumed error distribution and reordering them with inferred error ranks via the secant method, rather than estimating their values directly. The results of a synthetic case and a real case using the hydrologic models GR4J and HYMOD show several benefits of this new approach: 1) the reordering strategy via the secant method can significantly improve the efficiency and accuracy of input error quantification, and consequently promote the parameter estimation; 2) an autoregressive model can address the persistence of hydrologic residual errors in calibration; 3) the lag-time between the forcing data and the corresponding response can be explicitly acknowledged when the system exhibits a delayed response. This developed BEAR algorithm can be extended to other environmental modeling studies with correlated or/and delayed responses although its ability is limited by the impacts of model structural error and the output observational error.

Bayesian Frequency-Dependent AVO Inversion Using an Improved Markov Chain Monte Carlo Method for Quantitative Gas Saturation Prediction in a Thin Layer

One of the main objectives in the hydrocarbon reservoir characterization is determining rock and fluid properties that rely extensively on inference from seismic observations. In this letter, we present a novel Bayesian prestack inversion method using frequency-dependent amplitude versus offset (AVO) analysis with the goal to directly estimate gas saturation and porosity of a target thin reservoir zone. The proposed methodology is based on an improved Markov chain Monte Carlo (MCMC) sampling algorithm, which is computationally very coefficient due to its satisfactory acceptance probability and the convergence speed of Markov chains. Using a nonlinear rock physics model (RPM), properly calibrated for the investigating area, and a seismic forward operator based on the frequency-domain propagator matrix approach in the Bayesian inversion framework, we then evaluate the full posterior probability distribution of petrophysical parameters conditioned to seismic data and available prior information, using the MCMC algorithm in which we iteratively sample within the petrophysical property space. The proposed inversion approach is validated through applications to a synthetic reservoir model and the real seismic data from gas-bearing reservoirs with strong velocity dispersion.

Bayesian quantile semiparametric mixed-effects double regression models

Semiparametric mixed-effects double regression models have been used for analysis of longitudinal data in a variety of applications, as they allow researchers to jointly model the mean and variance of the mixed-effects as a function of predictors. However, these models are commonly estimated based on the normality assumption for the errors and the results may thus be sensitive to outliers and/or heavy-tailed data. Quantile regression is an ideal alternative to deal with these problems, as it is insensitive to heteroscedasticity and outliers and can make statistical analysis more robust. In this paper, we consider Bayesian quantile regression analysis for semiparametric mixed-effects double regression models based on the asymmetric Laplace distribution for the errors. We construct a Bayesian hierarchical model and then develop an efficient Markov chain Monte Carlo sampling algorithm to generate posterior samples from the full posterior distributions to conduct the posterior inference. The performance of the proposed procedure is evaluated through simulation studies and a real data application.

Heterogeneity learning for SIRS model: an application to the COVID-19

We propose a Bayesian Heterogeneity Learning approach for Susceptible-Infected-Removal-Susceptible (SIRS) model that allows underlying clustering patterns for transmission rate, recovery rate, and loss of immunity rate for the latest coronavirus (COVID-19) among different regions. Our proposed method provides simultaneously inference on parameter estimation and clustering information which contains both number of clusters and cluster configurations. Specifically, our key idea is to formulates the SIRS model into a hierarchical form and assign the Mixture of Finite mixtures priors for heterogeneity learning. The properties of the proposed models are examined and a Markov chain Monte Carlo sampling algorithm is used to sample from the posterior distribution. Extensive simulation studies are carried out to examine empirical performance of the proposed methods. We further apply the proposed methodology to analyze the state level COVID-19 data in U.S.

Incorporating compressor station multiple failure modes in risk evaluation of electricity-gas integrated energy systems

Multiple failure modes of gas compressor stations have not been modeled in the risk evaluation of electricity-gas integrated energy system (EGIES). This paper presents a method to incorporate the multiple failure modes of compressor stations with bypass system in the EGIES risk evaluation. Three outage models representing the multi-states and multi-failure-modes of compressor stations and a bi-level Monte Carlo sampling algorithm for the models are developed. A novel network model of EGIES considering the compressor stations with multiple failure modes is presented. An EGIES with a modified RBTS and an 11-node gas system is used to demonstrate the effectiveness of the proposed method and models. The results indicate that compressor stations have both positive and negative significant impacts on the risk of EGIES. The bypass systems for the compressor stations can effectively offset the negative impacts. Ignoring compressor stations and their station failures in the risk evaluation of EGIES may result in an inaccurate estimation of the system risk levels and even lead to a misleading conclusion in system planning.

S-wave velocity structures at Yedang Reservoir Dam inferred from amplification characteristics determined using H/V spectral ratios with background noise

Local shear wave (S-wave) velocity structure and sedimentary cover thickness are essential parameters that control local amplification of ground motion and associated seismic hazard during earthquakes. Recently, analysis of environmental background noise from individual stations has been used to estimate horizontal-to-vertical spectral ratio (H/V) curves. In the present study, we inverted H/V curves using environmental noise immediately prior to earthquake P-wave arrival to retrieve S-wave velocity (or shear velocity, V S) profiles at four temporary seismic monitoring stations (YDF, YDS, YDU, and YDD) near the Yedang Reservoir Dam. In the first step, we used a random search algorithm to constrain the subvolume of the parameter space (S-wave velocity structure) where the minimum of the misfit was located. In the second step, we independently applied two non-linear processes (a Monte Carlo sampling algorithm and a simulated annealing algorithm) to force the inversion towards an optimal solution, using the minimum misfit model determined in the first step as an initial estimate; we then compared the results. The feasibility and effectiveness of this two-step approach were verified by inversion of H/V curves for seismic noise recorded at four seismic monitoring stations near the Yedang Reservoir Dam. Borehole and topographical data from the four stations provided a well-constrained estimation of the local shear wave profile. Comparisons of synthetic and observed H/V curves showed that combining the two inversion algorithms efficiently overcame the extreme non-linearity of the inversion problem and provided a good resolution of S-wave structures at the Yedang Reservoir Dam. The S-wave velocity profile at the YDF station, which is situated on fresh, uniform bedrock, ranged from ∼2300 to 2700 m/s, which was consistent with the borehole data. Both the YDU station (which exhibited fundamental and first-order resonance frequency harmonics) and the YDD station (which exhibited fundamental, first-, and second-order resonance harmonics) showed significant stepwise velocity profiles for the entire depth. Thus, the entire depth can be regarded as three layers, including bedrock. At the YDS station, located near the dam spillway, the uppermost layer showed the lowest V S values among all four stations (∼726–728 m/s), suggesting that the concrete dam body is well weathered or partially infiltrated by water. Because of these low V S values of the concrete dam body, the Yedang Dam may require a full safety diagnosis in the near future. Analysis of the velocity profiles of all four stations near the Yedang Dam indicated that the bedrock, composed of fresh rock (∼V S > 2700 m/s), varies greatly in depth among stations.

Unbiased warped-area sampling for differentiable rendering

Differentiable rendering computes derivatives of the light transport equation with respect to arbitrary 3D scene parameters, and enables various applications in inverse rendering and machine learning. We present an unbiased and efficient differentiable rendering algorithm that does not require explicit boundary sampling. We apply the divergence theorem to the derivative of the rendering integral to convert the boundary integral into an area integral. We rewrite the converted area integral to a form that is suitable for Monte Carlo rendering. We then develop an efficient Monte Carlo sampling algorithm for solving the area integral. Our method can be easily plugged into a traditional path tracer and does not require dedicated data structures for sampling boundaries. We analyze the convergence properties through bias-variance metrics, and demonstrate our estimator's advantages over existing methods for some synthetic inverse rendering examples.

Compressed Monte Carlo with application in particle filtering

Bayesian models have become very popular over the last years in several fields such as signal processing, statistics, and machine learning. Bayesian inference requires the approximation of complicated integrals involving posterior distributions. For this purpose, Monte Carlo (MC) methods, such as Markov Chain Monte Carlo and importance sampling algorithms, are often employed. In this work, we introduce the theory and practice of a Compressed MC (C-MC) scheme to compress the statistical information contained in a set of random samples. In its basic version, C-MC is strictly related to the stratification technique, a well-known method used for variance reduction purposes. Deterministic C-MC schemes are also presented, which provide very good performance. The compression problem is strictly related to the moment matching approach applied in different filtering techniques, usually called as Gaussian quadrature rules or sigma-point methods. C-MC can be employed in a distributed Bayesian inference framework when cheap and fast communications with a central processor are required. Furthermore, C-MC is useful within particle filtering and adaptive IS algorithms, as shown by three novel schemes introduced in this work. Six numerical results confirm the benefits of the introduced schemes, outperforming the corresponding benchmark methods. A related code is also provided. (The code is provided athttp://www.lucamartino.altervista.org/CMC_CODE_pub_EX1.zip)

An Exact Auxiliary Variable Gibbs Sampler for a Class of Diffusions

Stochastic differential equations (SDEs) or diffusions are continuous-valued continuous-time stochastic processes widely used in the applied and mathematical sciences. Simulating paths from these processes is usually an intractable problem, and typically involves time-discretization approximations. We propose an exact Markov chain Monte Carlo sampling algorithm that involves no such time-discretization error. Our sampler is applicable to the problem of prior simulation from an SDE, posterior simulation conditioned on noisy observations, as well as parameter inference given noisy observations. Our work recasts an existing rejection sampling algorithm for a class of diffusions as a latent variable model, and then derives an auxiliary variable Gibbs sampling algorithm that targets the associated joint distribution. At a high level, the resulting algorithm involves two steps: simulating a random grid of times from an inhomogeneous Poisson process, and updating the SDE trajectory conditioned on this grid. Our work allows the vast literature of Monte Carlo sampling algorithms from the Gaussian process literature to be brought to bear to applications involving diffusions. We study our method on synthetic and real datasets, where we demonstrate superior performance over competing methods. Supplementary materials for this article are available online.

Joint modelling of longitudinal and survival data in the presence of competing risks with applications to prostate cancer data.

This research is motivated from the data from a large Selenium and Vitamin E Cancer Prevention Trial (SELECT). The prostate specific antigens (PSAs) were collected longitudinally, and the survival endpoint was the time to low-grade cancer or the time to high-grade cancer (competing risks). In this article, the goal is to model the longitudinal PSA data and the time-to-prostate cancer (PC) due to low- or high-grade. We consider the low-grade and high-grade as two competing causes of developing PC. A joint model for simultaneously analysing longitudinal and time-to-event data in the presence of multiple causes of failure (or competing risk) is proposed within the Bayesian framework. The proposed model allows for handling the missing causes of failure in the SELECT data and implementing an efficient Markov chain Monte Carlo sampling algorithm to sample from the posterior distribution via a novel reparameterization technique. Bayesian criteria, ΔDICSurv, and ΔWAICSurv, are introduced to quantify the gain in fit in the survival sub-model due to the inclusion of longitudinal data. A simulation study is conducted to examine the empirical performance of the posterior estimates as well as ΔDICSurv and ΔWAICSurv and a detailed analysis of the SELECT data is also carried out to further demonstrate the proposed methodology.

Bayesian hierarchical spatial regression models for spatial data in the presence of missing covariates with applications

AbstractIn many applications, survey data are collected from different survey centers in different regions. It happens that in some circumstances, response variables are completely observed while the covariates have missing values. In this article, we propose a joint spatial regression model for the response variable and missing covariates via a sequence of one‐dimensional conditional spatial regression models. We further construct a joint spatial model for missing covariate data mechanisms. The properties of the proposed models are examined and a Markov chain Monte Carlo sampling algorithm is used to sample from the posterior distribution. In addition, the Bayesian model comparison criteria, the modified deviance information criterion and the modified logarithm of the pseudo‐marginal likelihood, are developed to assess the fit of spatial regression models for spatial data. Extensive simulation studies are carried out to examine empirical performance of the proposed methods. We further apply the proposed methodology to analyze a real dataset from a Chinese Health and Nutrition Survey conducted in 2011.

Rao-Blackwellized sampling for batch and recursive Bayesian inference of Piecewise Affine models

This paper addresses batch (offline) and recursive (online) Bayesian inference of Piecewise Affine (PWA) regression models. By exploiting the particular structure of PWA models, efficient Rao-Blackwellized Monte Carlo sampling algorithms are developed to approximate the joint posterior distribution of the model parameters. Only the marginal posterior of the parameters used to describe the regressor-space partition is approximated, either in a batch mode using a Metropolis–Hastings Markov-Chain Monte Carlo (MCMC) sampler, or sequentially using particle filters, while the conditional distribution of the other model parameters is computed analytically. Probability distributions for the predicted outputs given new test inputs are derived and modifications of the proposed approaches to address maximum-a-posteriori estimate are discussed. The performance of the proposed algorithms is shown via a numerical example and through a benchmark case study on data-driven modelling of the electronic component placement process in a pick-and-place machine.

Uncertainty quantification for TRACE simulation of FIX-II No. 5052 test

The Best Estimate Plus Uncertainty approach requires the knowledge of input uncertainties for the uncertainty propagation with best-estimate codes. Inaccurate judgement of some model parameter uncertainties related to the dominant physical phenomena may result in misestimation of the safety margin. This paper presents a framework of inverse uncertainty quantification (UQ) to quantify model parameter uncertainties in order to address this issue. It is applied to TRACE simulation of a large break loss of coolant accident conducted on the FIX-II facility, and peak cladding temperature (PCT) is the simulation output. Sensitivity analysis identifies the parameters of the critical flow model as the most influential to the PCT. The inverse UQ is performed based on Bayesian framework, which adopts Markov Chain Monte Carlo sampling and surrogate modelling algorithms. The quantified uncertainties of the model parameters are the desired results from the inverse UQ process, which are useful in BEPU studies.

A Bayesian procedure for bandwidth selection in circular kernel density estimation

Abstract A Bayesian procedure for bandwidth selection in kernel circular density estimation is investigated, when the Markov chain Monte Carlo (MCMC) sampling algorithm is utilized for Bayes estimates. Under the quadratic and entropy loss functions, the proposed method is evaluated through a simulation study and real data sets, which were already discussed in the literature. The proposed Bayesian approach is very competitive in comparison with the existing classical global methods, namely plug-in and cross-validation techniques.

A Bayesian multi-risks survival (MRS) model in the presence of double censorings.

Semi-competing risks data include the time to a nonterminating event and the time to a terminating event, while competing risks data include the time to more than one terminating event. Our work is motivated by a prostate cancer study, which has one nonterminating event and two terminating events with both semi-competing risks and competing risks present as well as two censoring times. In this paper, we propose a new multi-risks survival (MRS) model for this type of data. In addition, the proposed MRS model can accommodate noninformative right-censoring times for nonterminating and terminating events. Properties of the proposed MRS model are examined in detail. Theoretical and empirical results show that the estimates of the cumulative incidence function for a nonterminating event may be biased if the information on a terminating event is ignored. A Markov chain Monte Carlo sampling algorithm is also developed. Our methodology is further assessed using simulations and also an analysis of the real data from a prostate cancer study. As a result, a prostate-specific antigen velocity greater than 2.0 ng/mL per year and higher biopsy Gleason scores are positively associated with a shorter time to death due to prostatecancer.

Second order linear differential equations with analytic uncertainties: Stochastic analysis via the computation of the probability density function

This paper concerns the analysis of random second order linear differential equations. Usually, solving these equations consists of computing the first statistics of the response process, and that task has been an essential goal in the literature. A more ambitious objective is the computation of the solution probability density function. We present advances on these two aspects in the case of general random non-autonomous second order linear differential equations with analytic data processes. The Fröbenius method is employed to obtain the stochastic solution in the form of a mean square convergent power series. We demonstrate that the convergence requires the boundedness of the random input coefficients. Further, the mean square error of the Fröbenius method is proved to decrease exponentially with the number of terms in the series, although not uniformly in time. Regarding the probability density function of the solution at a given time, which is the focus of the paper, we rely on the law of total probability to express it in closed-form as an expectation. For the computation of this expectation, a sequence of approximating density functions is constructed by reducing the dimensionality of the problem using the truncated power series of the fundamental set. We prove several theoretical results regarding the pointwise convergence of the sequence of density functions and the convergence in total variation. The pointwise convergence turns out to be exponential under a Lipschitz hypothesis. As the density functions are expressed in terms of expectations, we propose a crude Monte Carlo sampling algorithm for their estimation. This algorithm is implemented and applied on several numerical examples designed to illustrate the theoretical findings of the paper. After that, the efficiency of the algorithm is improved by employing the control variates method. Numerical examples corroborate the variance reduction of the Monte Carlo approach.