MCMC Sampling Research Articles

Borrowing strength and borrowing index for Bayesian hierarchical models

A novel borrowing strength measure and an overall borrowing index to characterize the strength of borrowing behaviors among subgroups are proposed for a given Bayesian hierarchical model. The constructions of the proposed indexes are based on the Mallow’s distance and can be easily computed using MCMC samples for univariate or multivariate posterior distributions. Consequently, the proposed indexes can serve as meaningful and useful exploratory tools to better understand the roles played by the priors in a hierarchical model, including their influences on the posteriors that are used to make statistical inferences. These relationships are otherwise ambiguous. The proposed methods can be applied to both the continuous and binary outcome variables. Furthermore, the proposed approach can be easily adapted to various settings of clinical trials, where Bayesian hierarchical models are deem appropriate. The effectiveness of the proposed method is illustrated using extensive simulation studies and a real data example.

Fast MCMC sampling algorithms on polytopes

We propose and analyze two new MCMC sampling algorithms, the Vaidya walk and the John walk, for generating samples from the uniform distribution over a polytope. Both random walks are sampling algorithms derived from interior point methods. The former is based on volumetric-logarithmic barrier introduced by Vaidya whereas the latter uses John's ellipsoids. We show that the Vaidya walk mixes in significantly fewer steps than the logarithmic-barrier based Dikin walk studied in past work. For a polytope in ℝd defined by n ≥ d linear constraints, we show that the mixing time from a warm start is bounded as O(n0.5d1.5), compared to the O(nd) mixing time bound for the Dikin walk. The cost of each step of the Vaidya walk is of the same order as the Dikin walk, and at most twice as large in terms of constant pre-factors. For the John walk, we prove an O(d2.5 log4(n/d) bound on its mixing time and conjecture that an improved variant of it could achieve a mixing time of O(d2 poly-log(n/d)). Additionally, we propose variants of the Vaidya and John walks that mix in polynomial time from a deterministic starting point. The speed-up of the Vaidya walk over the Dikin walk are illustrated in numerical examples.

ON STUDY OF EXPONENTIATED GAMMA DISTRIBUTION BASED ON UNIFIED HYBRID CENSORED DATA

In this article, we will study the estimation of the unknown parameters for exponentiated gamma distribution as well as a survival function, failure rate function and the coefficient of variation based on unified hybrid censored data. In addition that, we will study maximum likelihood and Bayesian estimates. To calculate the Bayes estimates of the model parameters will beused Markov chain Monte Carlo method (MCMC). Gibbs within the Metropolis-Hasting algorithm has been applied to generate MCMC samples from the posterior density function and calculate approximate confidence intervals for the unknown parameters, survival, failure rate functions and coefficient of variation. All resultsobtained are based on the balanced-squared error loss, balanced linear-exponential loss, and balanced general entropy loss functions. At the end of article, real data has been used to determine how the estimation method can be used in practice.

Is-ClusterMPP: clustering algorithm through point processes and influence space towards high-dimensional data

Clustering via marked point processes and influence space, Is-ClusterMPP, is a new unsupervised clustering algorithm through adaptive MCMC sampling of a marked point processes of interacting balls. The designed Gibbs energy cost function makes use of k-influence space information. It detects clusters of different shapes, sizes and unbalanced local densities. It aims at dealing also with high-dimensional datasets. By using the k-influence space, Is-ClusterMPP solves the problem of local heterogeneity in densities and prevents the impact of the global density in the detection of unbalanced classes. This concept reduces also the input values amount. The curse of dimensionality is handled by using a local subspace clustering principal embedded in a weighted similarity metric. Balls covering data points are constituting a configuration sampled from a marked point process (MPP). Due to the choice of the energy function, they tends to cover neighboring data, which share the same cluster. The statistical model of random balls is sampled through a Monte Carlo Markovian dynamical approach. The energy is balancing different goals. (1) The data driven objective function is provided according to k-influence space. Data in a high-dense region are favored to be covered by a ball. (2) An interaction part in the energy prevents the balls full overlap phenomenon and favors connected groups of balls. The algorithm through Markov dynamics, does converge towards configurations sampled from the MPP model. This algorithm has been applied in real benchmarks through gene expression data set of various sizes. Different experiments have been done to compare Is-ClusterMPP against the most well-known clustering algorithms and its efficiency is claimed.

Dynamic prediction using joint models of longitudinal and recurrent event data: a Bayesian perspective

In cardiovascular disease (CVD) studies, the events of interest may be recurrent (multiple occurrences from the same individual). During the study follow-up, longitudinal measurements are often available and these measurements are highly predictive of event recurrences. It is of great clinical interest to make personalized prediction of the next occurrence of recurrent events using the available clinical information, because it enables clinicians to make more informed and personalized decisions and recommendations. To this end, we propose a joint model of longitudinal and recurrent event data. We develop a Bayesian approach for model inference and a dynamic prediction framework for predicting target subjects' future outcome trajectories and risk of next recurrent event, based on their data up to the prediction time point. To improve computation efficiency, embarrassingly parallel MCMC (EP-MCMC) method is utilized. It partitions the data into multiple subsets, runs MCMC sampler on each subset, and applies random partition trees to combine the posterior draws from all subsets. Our method development is motivated by and applied to the Antihypertensive and Lipid-Lowering Treatment to Prevent Heart Attack Trial (ALLHAT), one of the largest CVD studies to compare the effectiveness of medications to treat hypertension.

Implicitly Localized MCMC Sampler to Cope With Non-local/Non-linear Data Constraints in Large-Size Inverse Problems

Many practical applications involve the resolution of large-size inverse problems, without providing more than a moderate-size sample to describe the prior probability distribution. In this situation, additional information must be supplied to augment the effective dimension of the available sample. This is the role played by covariance localization in the large-size applications of the ensemble Kalman filter. In this paper, it is suggested that covariance localization can also be efficiently applied to an approximate variant of the Metropolis/Hastings algorithm, by modulating the ensemble members by the large-scale patterns of other members. Modulation is used to design a (global) proposal probability distribution (i)~that can be sampled at a very low cost (proportional to the size of the state vector, with a small integer coefficient),(ii)~that automatically accounts for a localized prior covariance, and (iii)~that leads to an efficient sampler for the augmented prior probability distribution or for the posterior probability distribution. The resulting algorithm is applied to an academic example, illustrating (i)~the effectiveness of covariance localization, (ii)~the ability of the method to deal with nonlocal/nonlinear observation operators and non-Gaussian observation errors, (iii)~the possibility to deal with non-Gaussian (even discrete) prior marginal distributions, by including (stochastic) anamorphosis transformations, (iv)~the reliability, resolution and optimality of the updated ensemble, using probabilistic scores appropriate to a non-Gaussian posterior distribution, and (v)~the scalability of the algorithm as a function of the size of the problem. The evaluation of the computational complexity of the algorithm suggests that it could become numerically competitive with local ensemble Kalman filters, even in the absence of nonlocal constraints, especially if the localization radius is large. All codes necessary to reproduce the example application described in this paper are openly available from github.com/brankart/ensdam.

Euclid-era cosmology for everyone: Neural net assisted MCMC sampling for the joint 3x2 likelihood

Abstract We develop a fully non-invasive use of machine learning in order to enable open research on Euclid-sized data sets. Our algorithm leaves complete control over theory and data analysis, unlike many black-box like uses of machine learning. Focusing on a ‘3x2 analysis’ which combines cosmic shear, galaxy clustering and tangential shear at a Euclid-like sky coverage, we arrange a total of 348000 data points into data matrices whose structure permits not only an easy prediction by neural nets, but it additionally permits the essential removal from the data of patterns which the neural nets could not ‘understand’. The latter provides an often lacking mechanism to control and debias the inference of physics. The theoretical backbone to our neural net training can be any conventional (deterministic) theory code, where we chose CLASS. After training, we infer the seven parameters of a wCDM cosmology by Monte Carlo Markov sampling posteriors at Euclid-like precision within a day. We publicly provide the neural nets which memorise and output all 3x2 power spectra at a Euclid-like sky coverage and redshift binning.

A new method for determining the benchmark dose tolerable region and endpoint probabilities for toxicology experiments

ABSTRACT The increase of exposure to toxic materials and hazardous chemicals is a major concern due to the adverse effect on human health. Among the major concerns of toxicologists is to determine acceptable levels of exposure to hazardous substances. Current approaches often evaluate each endpoint and stressor individually. We propose a novel method to simultaneously determine the Benchmark Dose Tolerable Region (BMDTR) for multiple endpoint and multiple stressor studies by adopting a Bayesian approach. A main concern while assessing the combined toxicological effect of a chemical mixture is the anticipated type of the combined action (i.e. synergistic or antagonistic); thus it was essential to account for interaction effects to handle this situation, imposing more challenges due to the non-linearity of the tolerable region. The proposed method will be evaluated using two approaches, the first one using the estimated value of the posterior median and the second approach using all MCMC samples from the posterior distribution. Furthermore, we propose a new method to determine the endpoint probabilities for each endpoint, which reflects the importance of each endpoint in helping determining the boundaries of the benchmark dose tolerable region (BMDTR).

Estimation in the Complementary Exponential Geometric Distribution Based on Progressive Type-II Censored Data

Complementary exponential geometric distribution has many applications in survival and reliability analysis. Due to its importance, in this study, we are aiming to estimate the parameters of this model based on progressive type-II censored observations. To do this, we applied the stochastic expectation maximization method and Newton–Raphson techniques for obtaining the maximum likelihood estimates. We also considered the estimation based on Bayesian method using several approximate: MCMC samples, Lindely approximation and Metropolis–Hasting algorithm. In addition, we considered the shrinkage estimators based on Bayesian and maximum likelihood estimators. Then, the HPD intervals for the parameters are constructed based on the posterior samples from the Metropolis–Hasting algorithm. In the sequel, we obtained the performance of different estimators in terms of biases, estimated risks and Pitman closeness via Monte Carlo simulation study. This paper will be ended up with a real data set example for illustration of our purpose.

Fast likelihood-free cosmology with neural density estimators and active learning

Abstract Likelihood-free inference provides a framework for performing rigorous Bayesian inference using only forward simulations, properly accounting for all physical and observational effects that can be successfully included in the simulations. The key challenge for likelihood-free applications in cosmology, where simulation is typically expensive, is developing methods that can achieve high-fidelity posterior inference with as few simulations as possible. Density-estimation likelihood-free inference (DELFI) methods turn inference into a density estimation task on a set of simulated data-parameter pairs, and give orders of magnitude improvements over traditional Approximate Bayesian Computation approaches to likelihood-free inference. In this paper we use neural density estimators (NDEs) to learn the likelihood function from a set of simulated datasets, with active learning to adaptively acquire simulations in the most relevant regions of parameter space on-the-fly. We demonstrate the approach on a number of cosmological case studies, showing that for typical problems high-fidelity posterior inference can be achieved with just $\mathcal {O}(10^3)$ simulations or fewer. In addition to enabling efficient simulation-based inference, for simple problems where the form of the likelihood is known, DELFI offers a fast alternative to MCMC sampling, giving orders of magnitude speed-up in some cases. Finally, we introduce pydelfi – a flexible public implementation of DELFI with NDEs and active learning – available at https://github.com/justinalsing/pydelfi.

Revealing uninfected and infected target cell dynamics from peripheral blood data in highly and less pathogenic simian/human immunodeficiency virus infected Rhesus macaque

Since chimeric simian and human immunodeficiency viruses (SHIVs) used here, that is, SHIV-#64 and -KS661 utilize both CCR5 and CXCR4 chemokine receptors, they have broad target cell properties. A highly pathogenic SHIV strain, SHIV-KS661, causes an infection that systemically depletes the CD4+ T cells of Rhesus macaques, while a less pathogenic strain, SHIV-#64, does not cause severe symptoms in the macaques. In our previous studies, we established in vitro quantification system for virus infection dynamics, and concluded that SHIV-KS661 effectively produces infectious virions compared with SHIV-#64 in the HSC-F cell culture. However, in vivo dynamics of SHIV infection have not been well understood. To quantify SHIV-#64 and -KS661 infection dynamics in Rhesus macaques, we developed a novel approach and analyzed total CD4+ T cells and viral load in peripheral blood, and reproduced the expected dynamics for the uninfected and infected CD4+ T cells in silico. Using our previous cell culture experimental datasets, we revealed that an infection rate constant is different between SHIV-#64 and -KS661, but the viral production rate and the death rate are similar for the both strains. Thus, here, we assumed these relations in our in vivo data and carried out the data fitting. We performed Bayesian estimation for the whole dataset using MCMC sampling, and simultaneously fitted our novel model to total CD4+ T cells and viral load of SHIV-#64 and -KS661 infection. Our analyses explained that the Malthusian parameter (i.e., fitness of virus infection) and the basic reproduction number (i.e., potential of virus infection) for SHIV-KS661 are significantly higher than those of SHIV-#64. In addition, we demonstrated that the number of uninfected CD4+ T cells in SHIV-KS661 infected Rhesus macaques decreases to the significantly lower value than that before the inoculation several days earlier compared with SHIV-#64 infection. Taken together, the differences between SHIV-#64 and -KS661 infection before the peak viral load might determine the subsequent destiny, that is, whether the systemic CD4+ T cell depletion occurs or the host immune response develop.

Assessment of water quality in small agricultural ponds in Nara, Japan

In Japan, since ancient times, numerous small- to medium-size agricultural ponds have been constructed to irrigate paddy rice fields. Economic growth and urbanization in the modern era have resulted in a decreasing agricultural land area. With the aging of farmers in the country and due to decreasing management efforts, the degradation of the quality of water and the ponds’ environment has become an emerging issue in many regions. This study was conducted to investigate the status of the quality of water in agricultural ponds in the Yamato Plain, the center of agriculture, politics, economy, and culture of Nara Prefecture, Japan. Water quality data from 71 ponds were collected during 2013 to 2017. The water quality parameters of interest were pH, DO, EC, TN, TP, COD, SS, and Chl.a. Of these, pH, DO, and EC were measured at sampling sites using a portable apparatus, while water samples were collected for the analysis of TN, TP, COD, SS, and Chl.a in the laboratory. The results of the water quality analysis indicated deteriorated pond water for use of agriculture and also for the environment. Principal component analysis was applied to the water quality data for the assessment of the degradation level of the ponds. Results from the analysis showed a dominance of organic contamination. The MCMC sampling method was then applied to develop a Bayesian multiple regression model to quantify the range of correlation between the water quality indicator and pond catchment land use pattern. The model revealed pond water deteriorated with an increased residential area in the catchment, while an increased forest area tends to have a positive effect on pond water quality. The model has the potential to be a useful tool for the management of the environmental quality of ponds.

A tutorial on generalizing the default Bayesian t-test via posterior sampling and encompassing priors

With the advent of so-called default Bayesian hypothesis tests, scientists in applied fields have gained access to a powerful and principled method for testing hypotheses. However, such default tests usually come with a compromise, requiring the analyst to accept a one-size-fits-all approach to hypothesis testing. Further, such tests may not have the flexibility to test problems the scientist really cares about. In this tutorial, I demonstrate a flexible approach to generalizing one specific default test (the JZS t-test; Rouder et al., 2009) that is becoming increasingly popular in the social and behavioral sciences. The approach uses two theoretical results, the Savage-Dickey density ratio (Dickey and Lientz, 1980) and the technique of encompassing priors (Klugkist et al., 2005) in combination with MCMC sampling via an easy-to-use probabilistic modeling package for R called Greta. Through a comprehensive mathematical description of the techniques as well as illustrative examples, the reader is presented with a general, flexible workflow that can be extended to solve problems relevant to his or her own work.

Two Jovian Planets around the Giant Star HD 202696: A Growing Population of Packed Massive Planetary Pairs around Massive Stars?

Abstract We present evidence for a new two-planet system around the giant star HD 202696 (=HIP 105056, BD +26 4118). The discovery is based on public HIRES radial velocity (RV) measurements taken at Keck Observatory between 2007 July and 2014 September. We estimate a stellar mass of for HD 202696, which is located close to the base of the red giant branch. A two-planet self-consistent dynamical modeling MCMC scheme of the RV data followed by a long-term stability test suggests planetary orbital periods of P b = and P c = days, eccentricities of e b = and e c = , and minimum dynamical masses of m b = and m c = M Jup, respectively. Our stable MCMC samples are consistent with orbital configurations predominantly in a mean period ratio of 11:6 and its close-by high-order mean-motion commensurabilities with low eccentricities. For the majority of the stable configurations, we find an aligned or anti-aligned apsidal libration (i.e., Δω librating around 0° or 180°), suggesting that the HD 202696 system is likely dominated by secular perturbations near the high-order 11:6 mean-motion resonance. The HD 202696 system is yet another Jovian-mass pair around an intermediate-mass star with a period ratio below the 2:1 mean-motion resonance. Therefore, the HD 202696 system is an important discovery that may shed light on the primordial disk–planet properties needed for giant planets to break the strong 2:1 mean-motion resonance and settle in more compact orbits.

Building a telescope to look into high-dimensional image spaces

An image pattern can be represented by a probability distribution whose density is concentrated on different low-dimensional subspaces in the high-dimensional image space. Such probability densities have an astronomical number of local modes corresponding to typical pattern appearances. Related groups of modes can join to form macroscopic image basins that represent pattern concepts. Recent works use neural networks that capture high-order image statistics to learn Gibbs models capable of synthesizing realistic images of many patterns. However, characterizing a learned probability density to uncover the Hopfield memories of the model, encoded by the structure of the local modes, remains an open challenge. In this work, we present novel computational experiments that map and visualize the local mode structure of Gibbs densities. Efficient mapping requires identifying the global basins without enumerating the countless modes. Inspired by Grenander's jump-diffusion method, we propose a new MCMC tool called Attraction-Diffusion (AD) that can capture the macroscopic structure of highly non-convex densities by measuring metastability of local modes. AD involves altering the target density with a magnetization potential penalizing distance from a known mode and running an MCMC sample of the altered density to measure the stability of the initial chain state. Using a low-dimensional generator network to facilitate exploration, we map image spaces with up to 12,288 dimensions (64 $\times$ 64 pixels in RGB). Our work shows: (1) AD can efficiently map highly non-convex probability densities, (2) metastable regions of pattern probability densities contain coherent groups of images, and (3) the perceptibility of differences between training images influences the metastability of image basins.

MontePython 3: Boosted MCMC sampler and other features

Abstract MontePython is a parameter inference package for cosmology. We present the latest development of the code over the past couple of years. We explain, in particular, two new ingredients both contributing to improve the performance of Metropolis–Hastings sampling: an adaptation algorithm for the jumping factor, and a calculation of the inverse Fisher matrix, which can be used as a proposal density. We present several examples to show that these features speed up convergence and can save many hundreds of CPU-hours in the case of difficult runs, with a poor prior knowledge of the covariance matrix. We also summarize all the functionalities of MontePython in the current release, including new likelihoods and plotting options.

A rigorous EFT-based forward model for large-scale structure

Conventional approaches to cosmology inference from galaxy redshift surveys are based on n-point functions, which are under rigorous perturbative control on sufficiently large scales. Here, we present an alternative approach, which employs a likelihood at the level of the galaxy density field. By integrating out small-scale modes based on effective-field theory arguments, we prove that this likelihood is under perturbative control if certain specific conditions are met. We further show that the information captured by this likelihood is equivalent to the combination of the next-to-leading order galaxy power spectrum, leading-order bispectrum, and BAO reconstruction. Combined with MCMC sampling and MAP optimization techniques, our results allow for fully Bayesian cosmology inference from large-scale structure that is under perturbative control. We illustrate this via a first demonstration of unbiased cosmology inference from nonlinear large-scale structure using this likelihood. In particular, we show unbiased estimates of the power spectrum normalization σ8 from a catalog of simulated dark matter halos, where nonlinear information is crucial in breaking the b1−σ8 degeneracy.

Bayesian multiple logistic regression for case-control GWAS.

Genetic variants in genome-wide association studies (GWAS) are tested for disease association mostly using simple regression, one variant at a time. Standard approaches to improve power in detecting disease-associated SNPs use multiple regression with Bayesian variable selection in which a sparsity-enforcing prior on effect sizes is used to avoid overtraining and all effect sizes are integrated out for posterior inference. For binary traits, the logistic model has not yielded clear improvements over the linear model. For multi-SNP analysis, the logistic model required costly and technically challenging MCMC sampling to perform the integration. Here, we introduce the quasi-Laplace approximation to solve the integral and avoid MCMC sampling. We expect the logistic model to perform much better than multiple linear regression except when predicted disease risks are spread closely around 0.5, because only close to its inflection point can the logistic function be well approximated by a linear function. Indeed, in extensive benchmarks with simulated phenotypes and real genotypes, our Bayesian multiple LOgistic REgression method (B-LORE) showed considerable improvements (1) when regressing on many variants in multiple loci at heritabilities ≥ 0.4 and (2) for unbalanced case-control ratios. B-LORE also enables meta-analysis by approximating the likelihood functions of individual studies by multivariate normal distributions, using their means and covariance matrices as summary statistics. Our work should make sparse multiple logistic regression attractive also for other applications with binary target variables. B-LORE is freely available from: https://github.com/soedinglab/b-lore.

Simple Structure Detection Through Bayesian Exploratory Multidimensional IRT Models

In modern validity theory, a major concern is the construct validity of a test, which is commonly assessed through confirmatory or exploratory factor analysis. In the framework of Bayesian exploratory Multidimensional Item Response Theory (MIRT) models, we discuss two methods aimed at investigating the underlying structure of a test, in order to verify if the latent model adheres to a chosen simple factorial structure. This purpose is achieved without imposing hard constraints on the discrimination parameter matrix to address the rotational indeterminacy. The first approach prescribes a 2-step procedure. The parameter estimates are obtained through an unconstrained MCMC sampler. The simple structure is, then, inspected with a post-processing step based on the Consensus Simple Target Rotation technique. In the second approach, both rotational invariance and simple structure retrieval are addressed within the MCMC sampling scheme, by introducing a sparsity-inducing prior on the discrimination parameters. Through simulation as well as real-world studies, we demonstrate that the proposed methods are able to correctly infer the underlying sparse structure and to retrieve interpretable solutions.

PRISM revisited: Declarative implementation of a probabilistic programming language using multi-prompt delimited control

PRISM is a probabilistic programming language based on Prolog, augmented with primitives to represent probabilistic choice. It is implemented using a combination of low level support from a modified version of B-Prolog, source level program transformation, and libraries for inference and learning implemented in C. More recently, developers working with functional programming languages have taken the approach of embedding probabilistic primitives into an existing language, with little or no modification to the host language, often by using delimited continuations. Captured continuations represent pieces of the probabilistic program which can be manipulated to achieve a great variety of computational effects useful for inference. In this paper, I will describe an approach based on delimited control operators recently introduced into SWI Prolog. These are used to create a system of nested effect handlers which together implement a core functionality of PRISM—the building of explanation graphs—entirely in Prolog and using an order of magnitude less code. Other declarative programming tools, such as constraint logic programming, are used to implement tools for inference, such as the inside-outside and EM algorithms, lazy best-first explanation search, and MCMC samplers. By embedding the functionality of PRISM into SWI Prolog, users gain access to its rich libraries and development environment. By expressing the functionality of PRISM in a small amount of pure, high-level Prolog, this implementation facilitates further experimentation with the mechanisms of probabilistic logic programming, including new probabilistic modelling features and inference algorithms, such as variational inference in models with real-valued variables.