Markov Chain Monte Carlo Implementation Research Articles

The Orbit of Warm Jupiter WASP-106 b is Aligned with its Star

Understanding orbital obliquities, or the misalignment angles between a star’s rotation axis and the orbital axis of its planets, is crucial for unraveling the mechanisms of planetary formation and migration. In this study, we present an analysis of Rossiter–McLaughlin (RM) observations of the warm Jupiter exoplanet WASP-106 b. The high-precision radial velocity measurements were made with HARPS and HARPS-N during the transit of this planet. We aim to constrain the orientation of the planet’s orbit relative to its host star’s rotation axis. The RM observations are analyzed using a code which models the RM anomaly together with the Keplerian orbit given several parameters in combination with a Markov chain Monte Carlo implementation. We measure the projected stellar obliquity in the WASP-106 system for the first time and find λ = (−1 ± 11)°, supporting the theory of quiescent migration through the disk.

Layer adaptive node selection in Bayesian neural networks: Statistical guarantees and implementation details

Sparse deep neural networks have proven to be efficient for predictive model building in large-scale studies. Although several works have studied theoretical and numerical properties of sparse neural architectures, they have primarily focused on the edge selection. Sparsity through edge selection might be intuitively appealing; however, it does not necessarily reduce the structural complexity of a network. Instead pruning excessive nodes leads to a structurally sparse network with significant computational speedup during inference. To this end, we propose a Bayesian sparse solution using spike-and-slab Gaussian priors to allow for automatic node selection during training. The use of spike-and-slab prior alleviates the need of an ad-hoc thresholding rule for pruning. In addition, we adopt a variational Bayes approach to circumvent the computational challenges of traditional Markov Chain Monte Carlo (MCMC) implementation. In the context of node selection, we establish the fundamental result of variational posterior consistency together with the characterization of prior parameters. In contrast to the previous works, our theoretical development relaxes the assumptions of the equal number of nodes and uniform bounds on all network weights, thereby accommodating sparse networks with layer-dependent node structures or coefficient bounds. With a layer-wise characterization of prior inclusion probabilities, we discuss the optimal contraction rates of the variational posterior. We empirically demonstrate that our proposed approach outperforms the edge selection method in computational complexity with similar or better predictive performance. Our experimental evidence further substantiates that our theoretical work facilitates layer-wise optimal node recovery.

Bayesian modeling of the temporal evolution of seismicity using the ETAS.inlabru package

The epidemic type aftershock sequence (ETAS) model is widely used to model seismic sequences and underpins operational earthquake forecasting (OEF). However, it remains challenging to assess the reliability of inverted ETAS parameters for numerous reasons. For example, the most common algorithms just return point estimates with little quantification of uncertainty. At the same time, Bayesian Markov chain Monte Carlo implementations remain slow to run and do not scale well, and few have been extended to include spatial structure. This makes it difficult to explore the effects of stochastic uncertainty. Here, we present a new approach to ETAS modeling using an alternative Bayesian method, the integrated nested Laplace approximation (INLA). We have implemented this model in a new R-Package called ETAS.inlabru, which is built on the R packages R-INLA and inlabru. Our study has included extending these packages, which provided tools for modeling log-Gaussian Cox processes, to include the self-exciting Hawkes process that ETAS is a special case of. While we just present the temporal component here, the model scales to a spatio-temporal model and may include a variety of spatial covariates. This is a fast method that returns joint posteriors on the ETAS background and triggering parameters. Using a series of synthetic case studies, we explore the robustness of ETAS inversions using this method of inversion. We also included runnable notebooks to reproduce the figures in this article as part of the package's GitHub repository. We demonstrate that reliable estimates of the model parameters require that the catalog data contain periods of relative quiescence, as well as triggered sequences. We explore the robustness of the method under stochastic uncertainty in the training data and show that the method is robust to a wide range of starting conditions. We show how the inclusion of historic earthquakes prior to the modeled time window affects the quality of the inversion. Finally, we show that rate-dependent incompleteness of earthquake catalogs after large earthquakes have a significant and detrimental effect on the ETAS posteriors. We believe that the speed of the inlabru inversion, which includes a rigorous estimation of uncertainty, will enable a deeper exploration of how to use ETAS robustly for seismicity modeling and operational earthquake forecasting.

Laplace approximation for conditional autoregressive models for spatial data of diseases

Conditional autoregressive (CAR) distributions are used to account for spatial autocorrelation in small areal or lattice data to assess the spatial risks of diseases. The intrinsic CAR (ICAR) distribution has been primarily used as the priori distribution of spatially autocorrelated random variables in the framework of Bayesian statistics. The posterior distributions of spatial variates and unknown parameters of Bayesian ICAR models are estimated with the Markov chain Monte Carlo (MCMC) methods or integrated nested Laplace approximation (INLA), which may suffer from failures in numeric convergence. This study used the Laplace approximation, a fast computational method available in software Template Model Builder (TMB), for the maximum likelihood estimation (MLEs) of the ICAR model parameters. This study used the TMB to integrate out the latent spatial variates for the fast computations of marginal likelihood functions. This study compared the runtime and performance among the TMB, MCMC, and INLA implementations with three case studies of human diseases in the United Kingdom and the United States. The MLEs of the ICAR model with TMB were similar to those by the MCMC and INLA methods. The TMB implementation was faster than the MCMC (up to 100–200 times) and INLA (nine times) models.• This study built conditional autoregressive models in template model builder• TMB implementation was 100-200 times faster than the MCMC method• TMB implementation was also faster than Bayesian approximation with R INLA

Markov chain Monte Carlo for petrophysical inversion

Stochastic petrophysical inversion is a method used to predict reservoir properties from seismic data. Recent advances in stochastic optimization allow generating multiple realizations of rock and fluid properties conditioned on seismic data. To match the measured data and represent the uncertainty of the model variables, many realizations are generally required. Stochastic sampling and optimization of spatially correlated models are computationally demanding. Monte Carlo methods allow quantifying the uncertainty of the model variables but are impractical for high-dimensional models with spatially correlated variables. We have developed a Bayesian approach based on an efficient implementation of the Markov chain Monte Carlo (MCMC) method for the inversion of seismic data for the prediction of reservoir properties. Our Bayesian approach includes an explicit vertical correlation model in the proposal distribution. It is applied trace by trace, and the lateral continuity model is imposed by using the previously simulated values at the adjacent traces as conditioning data for simulating the initial model at the current trace. The methodology is first presented for a 1D problem to test the vertical correlation, and it is extended to 2D problems by including the lateral correlation and comparing two novel implementations based on sequential sampling. Our method is applied to synthetic data to estimate the posterior distribution of the petrophysical properties conditioned on the measured seismic data. The results are compared with an MCMC implementation without lateral correlation and demonstrate the advantage of integrating a spatial correlation model.

Quantile Regression Neural Networks: A Bayesian Approach

This article introduces a Bayesian neural network estimation method for quantile regression assuming an asymmetric Laplace distribution (ALD) for the response variable. It is shown that the posterior distribution for feedforward neural network quantile regression is asymptotically consistent under a misspecified ALD model. This consistency proof embeds the problem from density estimation domain and uses bounds on the bracketing entropy to derive the posterior consistency over Hellinger neighborhoods. This consistency result is shown in the setting where the number of hidden nodes grow with the sample size. The Bayesian implementation utilizes the normal-exponential mixture representation of the ALD density. The algorithm uses Markov chain Monte Carlo (MCMC) simulation technique - Gibbs sampling coupled with Metropolis-Hastings algorithm. We have addressed the issue of complexity associated with the afore-mentioned MCMC implementation in the context of chain convergence, choice of starting values, and step sizes. We have illustrated the proposed method with simulation studies and real data examples.

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.

Bayesian Analysis of Multifidelity Computer Models With Local Features and Nonnested Experimental Designs: Application to the WRF Model

ABSTRACT Motivated by a multi-fidelity Weather Research and Forecasting (WRF) climate model application where the available simulations are not generated based on hierarchically nested experimental design, we develop a new co-kriging procedure called augmented Bayesian treed co-kriging. The proposed procedure extends the scope of co-kriging in two major ways. We introduce a binary treed partition latent process in the multifidelity setting to account for nonstationary and potential discontinuities in the model outputs at different fidelity levels. Moreover, we introduce an efficient imputation mechanism which allows the practical implementation of co-kriging when the experimental design is nonhierarchically nested by enabling the specification of semiconjugate priors. Our imputation strategy allows the design of an efficient reversible jump Markov chain Monte Carlo implementation that involves collapsed blocks and direct simulation from conditional distributions. We develop the Monte Carlo recursive emulator which provides a Monte Carlo proxy for the full predictive distribution of the model output at each fidelity level, in a computationally feasible manner. The performance of our method is demonstrated on benchmark examples and used for the analysis of a large-scale climate modeling application which involves the WRF model.

Importance sampling type estimators based on approximate marginal Markov chain Monte Carlo

AbstractWe consider importance sampling (IS) type weighted estimators based on Markov chain Monte Carlo (MCMC) targeting an approximate marginal of the target distribution. In the context of Bayesian latent variable models, the MCMC typically operates on the hyperparameters, and the subsequent weighting may be based on IS or sequential Monte Carlo (SMC), but allows for multilevel techniques as well. The IS approach provides a natural alternative to delayed acceptance (DA) pseudo‐marginal/particle MCMC, and has many advantages over DA, including a straightforward parallelization and additional flexibility in MCMC implementation. We detail minimal conditions which ensure strong consistency of the suggested estimators, and provide central limit theorems with expressions for asymptotic variances. We demonstrate how our method can make use of SMC in the state space models context, using Laplace approximations and time‐discretized diffusions. Our experimental results are promising and show that the IS‐type approach can provide substantial gains relative to an analogous DA scheme, and is often competitive even without parallelization.

A reduced order approach for probabilistic inversions of 3-D magnetotelluric data I: general formulation

SUMMARYSimulation-based probabilistic inversions of 3-D magnetotelluric (MT) data are arguably the best option to deal with the nonlinearity and non-uniqueness of the MT problem. However, the computational cost associated with the modelling of 3-D MT data has so far precluded the community from adopting and/or pursuing full probabilistic inversions of large MT data sets. In this contribution, we present a novel and general inversion framework, driven by Markov Chain Monte Carlo (MCMC) algorithms, which combines (i) an efficient parallel-in-parallel structure to solve the 3-D forward problem, (ii) a reduced order technique to create fast and accurate surrogate models of the forward problem and (iii) adaptive strategies for both the MCMC algorithm and the surrogate model. In particular, and contrary to traditional implementations, the adaptation of the surrogate is integrated into the MCMC inversion. This circumvents the need of costly offline stages to build the surrogate and further increases the overall efficiency of the method. We demonstrate the feasibility and performance of our approach to invert for large-scale conductivity structures with two numerical examples using different parametrizations and dimensionalities. In both cases, we report staggering gains in computational efficiency compared to traditional MCMC implementations. Our method finally removes the main bottleneck of probabilistic inversions of 3-D MT data and opens up new opportunities for both stand-alone MT inversions and multi-observable joint inversions for the physical state of the Earth’s interior.

Two-group Poisson-Dirichlet mixtures for multiple testing.

The simultaneous testing of multiple hypotheses is common to the analysis of high-dimensional data sets. The two-group model, first proposed by Efron, identifies significant comparisons by allocating observations to a mixture of an empirical null and an alternative distribution. In the Bayesian nonparametrics literature, many approaches have suggested using mixtures of Dirichlet Processes in the two-group model framework. Here, we investigate employing mixtures of two-parameter Poisson-Dirichlet Processes instead, and show how they provide a more flexible and effective tool for large-scale hypothesis testing. Our model further employs nonlocal prior densities to allow separation between the two mixture components. We obtain a closed-form expression for the exchangeable partition probability function of the two-group model, which leads to a straightforward Markov Chain Monte Carlo implementation. We compare the performance of our method for large-scale inference in a simulation study and illustrate its use on both a prostate cancer data set and a case-control microbiome study of the gastrointestinal tracts in children from underdeveloped countries who have been recently diagnosed with moderate-to-severediarrhea.

Households Health Status in Saudi Arabia: Spatial Distribution and Association With Socioeconomic Factors

This paper examines the spatial relationship between Saudi and non-Saudi people's health status and the socioeconomic composition of the neighbourhoods in which they live. Data were recorded from the National Population Health Survey (NPHS) performed by the Saudi General Authority for Statistics (GAS) in 2018. The survey counts 23,980,846 inhabitants grouped into 24,012 households who assessed their health status by gender and administrative region. Only people who are fifteen years of age and over and claiming poor health status were retained in the analysis. We used a Generalized Linear Spatial Model (GLSM) to study the relationship between Saudi and non-Saudi household’s health status and socioeconomic factors. A Gaussian process with a powered exponential spatial correlation function was introduced on the right-hand side of the model to consider the unexplained spatial variation in the data. The statistical results show the progressive increase in the number of Saudi and non-Saudi households claiming poor health status with the high Saudi unemployment rate, low average monthly income and high current daily smokers. The results of the statistical analyses show the wider potential of GLSM for analyzing data of this kind and the important risk of misleading interpretations when the non-spatial analysis is used on spatially structured data. The method of inference was Bayesian using Markov Chain Monte Carlo Implementation.

On the Anatomy of MCMC-Based Maximum Likelihood Learning of Energy-Based Models

This study investigates the effects of Markov chain Monte Carlo (MCMC) sampling in unsupervised Maximum Likelihood (ML) learning. Our attention is restricted to the family of unnormalized probability densities for which the negative log density (or energy function) is a ConvNet. We find that many of the techniques used to stabilize training in previous studies are not necessary. ML learning with a ConvNet potential requires only a few hyper-parameters and no regularization. Using this minimal framework, we identify a variety of ML learning outcomes that depend solely on the implementation of MCMC sampling.On one hand, we show that it is easy to train an energy-based model which can sample realistic images with short-run Langevin. ML can be effective and stable even when MCMC samples have much higher energy than true steady-state samples throughout training. Based on this insight, we introduce an ML method with purely noise-initialized MCMC, high-quality short-run synthesis, and the same budget as ML with informative MCMC initialization such as CD or PCD. Unlike previous models, our energy model can obtain realistic high-diversity samples from a noise signal after training.On the other hand, ConvNet potentials learned with non-convergent MCMC do not have a valid steady-state and cannot be considered approximate unnormalized densities of the training data because long-run MCMC samples differ greatly from observed images. We show that it is much harder to train a ConvNet potential to learn a steady-state over realistic images. To our knowledge, long-run MCMC samples of all previous models lose the realism of short-run samples. With correct tuning of Langevin noise, we train the first ConvNet potentials for which long-run and steady-state MCMC samples are realistic images.

On the use of approximate Bayesian computation Markov chain Monte Carlo with inflated tolerance and post-correction

SummaryApproximate Bayesian computation enables inference for complicated probabilistic models with intractable likelihoods using model simulations. The Markov chain Monte Carlo implementation of approximate Bayesian computation is often sensitive to the tolerance parameter: low tolerance leads to poor mixing and large tolerance entails excess bias. We propose an approach that involves using a relatively large tolerance for the Markov chain Monte Carlo sampler to ensure sufficient mixing and post-processing the output, leading to estimators for a range of finer tolerances. We introduce an approximate confidence interval for the related post-corrected estimators and propose an adaptive approximate Bayesian computation Markov chain Monte Carlo algorithm, which finds a balanced tolerance level automatically based on acceptance rate optimization. Our experiments show that post-processing-based estimators can perform better than direct Markov chain Monte Carlo targeting a fine tolerance, that our confidence intervals are reliable, and that our adaptive algorithm leads to reliable inference with little user specification.

Bayesian approach for proportional hazards mixture cure model allowing non-curable competing risk

ABSTRACTWith advancements in medical research, more and more diseases may be curable, which indicates some patients may not die of the disease of interest. Mixture cure models, which can capture patients being cured, attracts an increasing attention in practice. However, the existing mixture cure models only focused on the major event with a potential cure while ignoring the potential risk from other non-curable competing events, which are commonly seen in the real world. In this paper, we develop a Bayesian approach to estimate a proportional hazards mixture cure (PHMC) model allowing non-curable competing risk. Data augmentation method with latent binary cure indicators and event indicators are adopted to simplify the Markov chain Monte Carlo implementation. The baseline cumulative hazards for the PHMC model are formulated by counting processes with gamma process priors. Its performance is demonstrated through comprehensive simulation studies. Finally, the proposed method is applied to a prostate cancer clinical trial data set.

Comparison of Markov Chain Monte Carlo Software for the Evolutionary Analysis of Y-Chromosomal Microsatellite Data

The evolutionary analysis of genetic data is an important subject of modern bioscience, with practical applications in diverse fields. Parameters of interest in this context include effective population sizes, mutation rates, population growth rates and the times to most recent common ancestors. Studying Y-chromosomal microsatellite data, in particular, has proven useful to unravel the recent patrilineal history of Homo sapiens populations. We compared the individual analysis options and technical details of four software tools that are widely used for this purpose, namely BATWING, BEAST, IMa2 and LAMARC, all of which use Bayesian coalescent-based Markov chain Monte Carlo (MCMC) methods for parameter estimation. More specifically, we simulated datasets for either eight or 20 hypothetical Y-chromosomal microsatellites, assuming a mutation rate of 0.0030 per generation and a constant or exponentially increasing population size, and used these data to evaluate the parameter estimation capacity of each tool. The datasets comprised between 100 and 1000 samples. In addition to runtime, the practical utility of the tools of interest can also be expected to depend critically upon the convergence behavior of the actual MCMC implementation. In fact, we found that runtime increased, and convergence rate decreased, with increasing sample size as expected. BATWING performed best with respect to runtime and convergence behavior, but only supports simple evolutionary models. As regards the spectrum of evolutionary models covered, and also in terms of cross-platform usability, BEAST provided the greatest flexibility. Finally, IMa2 and LAMARC turned out best to incorporate elaborate migration models in the analysis process.

Multimodal Layered Transdimensional Inversion of Seismic Dispersion Curves With Depth Constraints

AbstractMuLTI (Multimodal Layered Transdimensional Inversion) is a Markov chain Monte Carlo implementation of Bayesian inversion for the probability distribution of shear wave velocity (Vs) as a function of depth. Based on Multichannel Analysis of Surface Wave methods, it requires as data (i) a Rayleigh‐wave dispersion curve and (ii) additional layer depth constraints, the latter we show significantly improve resolution compared to conventional unconstrained inversions. Such depth constraints may be drawn from any source (e.g., boreholes, complementary geophysical data) provided they also represent a seismic interface. We apply MuLTI to a Norwegian glacier, Midtdalsbreen, in which ground‐penetrating radar was used to constrain internal layers of snow, ice, and subglacial sediments, with transitions at 2 and 25.5 m, and whose Vs is assumed to be in the range 500–1,700, 1,700–1,950, and 200–2,800 m/s, respectively. Synthetic modeling demonstrates that MuLTI recovers the true model of Vs variation with depth. Significantly, compared to inversions without depth constraints, in this synthetic case MuLTI not only reduces by about a factor of 10 the error between the true and the best fitting model, but also reduces by a factor of 2 the vertically averaged spread of the distribution of Vs based on the 95% credible intervals. We further show that using frequencies above about 100 Hz lead to unconverged solutions due to mode ambiguities associated with fine spatial structures. For our acquired data on Midtdalsbreen, we use 14‐100 Hz data for which MuLTI produces a large‐scale converged inversion.

MODEL REGRESI LINIER BAYESIAN DENGAN APLIKASI PADA DATA PENUNDAAN PENERBANGAN

Bayesian linear regression is an approach to linear regression where statistical analysis depend of Bayesian inference. The Bayesian model on big data uses a summary of data statistics as input; Statistical summary can be calculated from each subset, then a statistical summary of the full dataset is obtained from the sum of the summary statistics for each subset. Recent developments in data science and research, produce large datasets that are too large to be analyzed as a whole due to the limitations of computer memory or storage capacity. To overcome this, a program package was introduced from R namely BayesSummaryStatLM for the Bayesian linear regression model with the Markov Chain Monte Carlo implementation that overcomes this limitation. Then the program package from R, ff is used to read data in large datasets while calculating statistics summary. In this study Bayesian linear regression model used with several choices of prior distribution for unknown model parameters, and illustrates in simulation data and real datasets for flight delay data in US 2008. The application of simulation data and flight delay data produces a plot of density functions for the β parameters has a shape resembling a plot of Normal distribution density function, whereas for plot parameters the density function has a shape resembling the plot of Inverse Gamma distribution density function. In the simulation data, the estimator for each parameter produced has a value that approach to the value of the specified parameter (True Value). This is also indicated by the narrow credible interval for each parameters.

Response to Cho and Liu, “Sampling from complicated and unknown distributions: Monte Carlo and Markov chain Monte Carlo methods for redistricting”

A question of legal significance is whether an enacted map of political districts is “typical.” Recent work has used Markov chain Monte Carlo (MCMC) methods to produce null distributions of maps in order to answer this question. A recent article by Cho and Liu critiques one particular implementation of MCMC for redistricting, that of Fifield et al. The goal of the present commentary is to draw attention to two facts omitted by Cho and Liu that, if included, would have severely weakened their conclusions. In particular, Cho and Liu point out that Fifield et al.’s algorithm fails to approximate a known target distribution, but neglect Fifield et al.’s use of parallel and simulated tempering, which greatly improves the approximation. Secondly, Cho and Liu argue that it is overly difficult to detect when Markov chains have mixed; they neglect to mention diagnostics used for this exact purpose in Fifield et al.

A Bayesian space-time model for clustering areal units based on their disease trends.

SummaryPopulation-level disease risk across a set of non-overlapping areal units varies in space and time, and a large research literature has developed methodology for identifying clusters of areal units exhibiting elevated risks. However, almost no research has extended the clustering paradigm to identify groups of areal units exhibiting similar temporal disease trends. We present a novel Bayesian hierarchical mixture model for achieving this goal, with inference based on a Metropolis-coupled Markov chain Monte Carlo ((MC)n}{}^3) algorithm. The effectiveness of the (MC)n}{}^3 algorithm compared to a standard Markov chain Monte Carlo implementation is demonstrated in a simulation study, and the methodology is motivated by two important case studies in the United Kingdom. The first concerns the impact on measles susceptibility of the discredited paper linking the measles, mumps, and rubella vaccination to an increased risk of Autism and investigates whether all areas in the Scotland were equally affected. The second concerns respiratory hospitalizations and investigates over a 10 year period which parts of Glasgow have shown increased, decreased, and no change in risk.