Unbiasing Hamiltonian Monte Carlo Algorithms for a General Hamiltonian Function

The Hamiltonian Monte Carlo (HMC) algorithm is a Markov Chain Monte Carlo (MCMC) technique, which combines the advantages of Hamiltonian dynamics methods and Metropolis Monte Carlo approach, to sample from complex distributions. The HMC algorithm incorporates gradient information in the dynamic trajectories and thus suppresses the random walk nature in traditional Markov chain simulation methods. This ensures rapid mixing, faster convergence, and improved efficiency of the Markov chain. The leapfrog method is generally used in discrete simulation of the dynamic transitions. In this paper, we refer to this as the leapfrog–HMC. The primary goal of this paper is to present the HMC algorithm as a tool for rapid sampling of high dimensional and complex distributions, and demonstrate its advantages over the classical Metropolis Monte Carlo technique. We demonstrate that the use of an adaptive–step discretization scheme in simulating the dynamic transitions results in an algorithm which significantly outperforms the leapfrog–HMC algorithm. Relevance to reservoir parameter estimation and uncertainty quantification will be discussed.

A reliable assessment of the posterior uncertainties is a crucial aspect of any amplitude versus angle (AVA) inversion due to the severe ill-conditioning of this inverse problem. To accomplish this task, numerical Markov chain Monte Carlo algorithms are usually used when the forward operator is nonlinear. The downside of these algorithms is the considerable number of samples needed to attain stable posterior estimations especially in high-dimensional spaces. To overcome this issue, we assessed the suitability of Hamiltonian Monte Carlo (HMC) algorithm for nonlinear target- and interval-oriented AVA inversions for the estimation of elastic properties and associated uncertainties from prestack seismic data. The target-oriented approach inverts the AVA responses of the target reflection by adopting the nonlinear Zoeppritz equations, whereas the interval-oriented method inverts the seismic amplitudes along a time interval using a 1D convolutional forward model still based on the Zoeppritz equations. HMC uses an artificial Hamiltonian system in which a model is viewed as a particle moving along a trajectory in an extended space. In this context, the inclusion of the derivative information of the misfit function makes possible long-distance moves with a high probability of acceptance from the current position toward a new independent model. In our application, we adopt a simple Gaussian a priori distribution that allows for an analytical inclusion of geostatistical constraints into the inversion framework, and we also develop a strategy that replaces the numerical computation of the Jacobian with a matrix operator analytically derived from a linearization of the Zoeppritz equations. Synthetic and field data inversions demonstrate that the HMC is a very promising approach for Bayesian AVA inversion that guarantees an efficient sampling of the model space and retrieves reliable estimations and accurate uncertainty quantifications with an affordable computational cost.

We show how the Hamiltonian Monte Carlo algorithm can sometimes be speeded up by “splitting” the Hamiltonian in a way that allows much of the movement around the state space to be done at low computational cost. One context where this is possible is when the log density of the distribution of interest (the potential energy function) can be written as the log of a Gaussian density, which is a quadratic function, plus a slowly-varying function. Hamiltonian dynamics for quadratic energy functions can be analytically solved. With the splitting technique, only the slowly-varying part of the energy needs to be handled numerically, and this can be done with a larger stepsize (and hence fewer steps) than would be necessary with a direct simulation of the dynamics. Another context where splitting helps is when the most important terms of the potential energy function and its gradient can be evaluated quickly, with only a slowly-varying part requiring costly computations. With splitting, the quick portion can be handled with a small stepsize, while the costly portion uses a larger stepsize. We show that both of these splitting approaches can reduce the computational cost of sampling from the posterior distribution for a logistic regression model, using either a Gaussian approximation centered on the posterior mode, or a Hamiltonian split into a term that depends on only a small number of critical cases, and another term that involves the larger number of cases whose influence on the posterior distribution is small.

The subordinated Poisson process is a generalization of the Poisson process, where the time of the original process is modelled using a stochastic process with non-decreasing paths. Using latent variables, the subordination provides additional flexibility, enabling modelling count data with larger jumps. The article discusses the parameter estimation of the subordinated Poisson process with Gamma and heavy-tailed Inverse Gaussian subordinators in Bayesian framework. The identifiability issues in the model parameters are analytically derived, and the models are then re-parameterized for identifiability. The methodology to infer the parameters for the subordinated Poisson process of order k is also proposed. Computation, using the Hamiltonian Monte Carlo algorithm, is demonstrated with detailed simulation studies. The accuracy of the inference procedure is established based on the low bias in the estimates and good convergence diagnostics. Data analyses are used to illustrate the proposed method.

We propose a novel approach to the detection of point-like sources of high-energy neutrinos. Motivated by evidence for emerging sources in existing data, we focus on the characterization and interpretation of these sources rather than the rejection of the background-only hypothesis. The hierarchical Bayesian model is implemented in the Stan platform, enabling computation of the posterior distribution with a Hamiltonian Monte Carlo algorithm. We simulate a population of weak neutrino sources detected by the IceCube experiment and use the resulting data set to demonstrate and validate our framework. We show that even for the challenging case of sources at the threshold of detection and using limited prior information, it is possible to correctly infer the source properties. Additionally, we demonstrate how modeling flexible connections between similar sources can be used to recover the contribution of sources that would not be detectable individually. While a direct comparison of our method to existing approaches is challenged by the fundamental differences in frequentist and Bayesian frameworks, we draw parallels where possible. In particular, we highlight how including more complexity into the source modeling can increase the sensitivity to sources and their populations.

ABSTRACTThe marginal regression model is convenient for analyzing correlated responses, including repeated measures and longitudinal data. This paper proposes a robust marginal linear model for analyzing a vector of univariate responses with correlated components by incorporating an innovative multivariate Huber distribution. It employs a flexible parameterization using modified Cholesky decomposition, provides a convenient approach for estimating the covariance matrix, and allows for subject‐varying the tuning parameter. Our research introduces a method for estimating parameters by employing the exact likelihood function through the Hamiltonian Monte Carlo algorithm. To highlight the advantage of our model, we carried out a simulation experiment and reanalyzed two real‐world case studies in the health and economics fields. The results indicate that our model offers a more robust analysis by assigning appropriate weights to extreme observations, thereby handling outliers more effectively than traditional models.

The standard errors (SE) of the maximum likelihood estimates (MLE) of the population parameter vector in nonlinear mixed effect models (NLMEM) are usually estimated using the inverse of the Fisher information matrix (FIM). However, at a finite distance, i.e. far from the asymptotic, the FIM can underestimate the SE of NLMEM parameters. Alternatively, the standard deviation of the posterior distribution, obtained in Stan via the Hamiltonian Monte Carlo algorithm, has been shown to be a proxy for the SE, since, under some regularity conditions on the prior, the limiting distributions of the MLE and of the maximum a posterior estimator in a Bayesian framework are equivalent. In this work, we develop a similar method using the Metropolis-Hastings (MH) algorithm in parallel to the stochastic approximation expectation maximisation (SAEM) algorithm, implemented in the saemix R package. We assess this method on different simulation scenarios and data from a real case study, comparing it to other SE computation methods. The simulation study shows that our method improves the results obtained with frequentist methods at finite distance. However, it performed poorly in a scenario with the high variability and correlations observed in the real case study, stressing the need for calibration.

Copula-based reliability analysis of degrading systems with dependent failures

In statistical learning theory, generalization error is used to quantify the degree to which a supervised machine learning algorithm may overfit to training data. Recent work [Xu and Raginsky (2017)] has established a bound on the generalization error of empirical risk minimization based on the mutual information $I$ ( $S$ ; W) between the algorithm input $S$ and the algorithm output W, when the loss function is sub-Gaussian. We leverage these results to derive generalization error bounds for a broad class of iterative algorithms that are characterized by bounded, noisy updates with Markovian structure. Our bounds are very general and are applicable to numerous settings of interest, including stochastic gradient Langevin dynamics (SGLD) and variants of the stochastic gradient Hamiltonian Monte Carlo (SGHMC) algorithm. Furthermore, our error bounds hold for any output function computed over the path of iterates, including the last iterate of the algorithm or the average of subsets of iterates, and also allow for non-uniform sampling of data in successive updates of the algorithm.

<strong class="journal-contentHeaderColor">Abstract.</strong> We present an efficient probabilistic workflow for the estimation of source parameters of induced seismic events in three-dimensional heterogeneous media. Our workflow exploits a linearized variant of the Hamiltonian Monte Carlo (HMC) algorithm. Compared to traditional Markov-Chain Monte Carlo (MCMC) algorithms, HMC is highly efficient in sampling high-dimensional model spaces. Through a linearization of the forward problem around the prior mean (i.e., the "best" initial model), this efficiency can be further improved. We show, however, that this linearization leads to a performance in which the output of an HMC chain strongly depends on the quality of the prior; in particular, because not all (induced) earthquake model parameters have a linear relationship with the recordings observed at the surface. To mitigate the importance of an accurate prior, we integrate the linearized HMC scheme into a workflow that (i) allows for a weak prior through linearization around various (initial) centroid locations, (ii) is able to converge to the mode containing the model with the (global) minimum misfit by means of an iterative HMC approach, and (iii) uses variance reduction as a criterion to include the output of individual Markov chains in the estimation of the posterior probability. Using a three-dimensional heterogeneous subsurface model of the Groningen gas field, we simulate an induced earthquake to test our workflow. We then demonstrate the virtue of our workflow by estimating the event's centroid (three parameters), moment tensor (six parameters), and the earthquake's origin time. We find that our workflow is able to recover the posterior probability of these source parameters rather well, even when the prior model information is inaccurate, imprecise, or both inaccurate and imprecise.

Models of Stellar Population Synthesis (SPS) provide a predictive framework for the spectral energy distribution (SED) of a galaxy. SPS predictions can be computationally intensive, creating a bottleneck for attempts to infer the physical properties of large populations of individual galaxies from their SEDs and photometry; these computational challenges are especially daunting for near-future cosmology surveys that will measure the photometry of billions of galaxies. In this paper, we explore a range of computational techniques aimed at accelerating SPS predictions of galaxy photometry using the JAX library to target GPUs. We study a particularly advantageous approximation to the calculation of galaxy photometry that speeds up the computation by a factor of 50 relative to the exact calculation. We introduce a novel technique for incorporating burstiness into models of galaxy star formation history that captures very short-timescale fluctuations with negligible increase in computation time. We study the performance of Hamiltonian Monte Carlo (HMC) algorithms in which individual chains are parallelized across independent GPU threads, finding that our pipeline can carry out Bayesian inference at a rate of approximately 1,000 galaxy posteriors per minute on a single GPU. Our results provide an update to standard benchmarks in the literature on the computational demands of SPS inference; our publicly available code enables previously-impractical Bayesian analyses of large galaxy samples, and includes several standalone modules that could be adopted to speedup existing SPS pipelines.

Importance sampling (IS) is a powerful Monte Carlo (MC) methodology for approximating integrals, for instance in the context of Bayesian inference. In IS, the samples are simulated from the so-called proposal distribution, and the choice of this proposal is key for achieving a high performance. In adaptive IS (AIS) methods, a set of proposals is iteratively improved. AIS is a relevant and timely methodology although many limitations remain yet to be overcome, e.g., the curse of dimensionality in high-dimensional and multi-modal problems. Moreover, the Hamiltonian Monte Carlo (HMC) algorithm has become increasingly popular in machine learning and statistics. HMC has several appealing features such as its exploratory behavior, especially in high-dimensional targets, when other methods suffer. In this letter, we introduce the novel Hamiltonian adaptive importance sampling (HAIS) method. HAIS implements a two-step adaptive process with parallel HMC chains that cooperate at each iteration. The proposed HAIS efficiently adapts a population of proposals, extracting the advantages of HMC. HAIS can be understood as a particular instance of the generic layered AIS family with an additional resampling step. HAIS achieves a significant performance improvement in high-dimensional problems w.r.t. state-of-the-art algorithms. We discuss the statistical properties of HAIS and show its high performance in two challenging examples.

Markov chain Monte Carlo (MCMC) algorithms have become powerful tools for Bayesian inference. However, they do not scale well to large-data problems. Divide-and-conquer strategies, which split the data into batches and, for each batch, run independent MCMC algorithms targeting the corresponding subposterior, can spread the computational burden across a number of separate computer cores. The challenge with such strategies is in recombining the subposteriors to approximate the full posterior. By creating a Gaussian-process approximation for each log-subposterior density we create a tractable approximation for the full posterior. This approximation is exploited through three methodologies: firstly a Hamiltonian Monte Carlo algorithm targeting the expectation of the posterior density provides a sample from an approximation to the posterior; secondly, evaluating the true posterior at the sampled points leads to an importance sampler that, asymptotically, targets the true posterior expectations; finally, an alternative importance sampler uses the full Gaussian-process distribution of the approximation to the log-posterior density to re-weight any initial sample and provide both an estimate of the posterior expectation and a measure of the uncertainty in it.

BackgroundThe goal of this study was to apply a dynamic modeling approach to assess the impact of COVID-19 control measures in nursing homes. Utah's strategy for COVID control in nursing homes included serial point prevalence surveys, triggered when a staff or resident case was detected. However, implementation of the policy varied across nursing homes.MethodsWe used state-collected testing data from 70 different Utah nursing homes from March 2020 to April 2021 to estimate a facility specific transmission. Data included facility id, COVID-19 test results, number of days since the first covid test administration, and a December 2020 population census. We constructed a stochastic Reed Frost based model, using the stan library in R using the Hamiltonian Monte Carlo algorithm to estimate parameter values.ResultsCOVID attack rates tended to peak between 50-100 days from the onset of included data, averaging 14.2% (Figure 1). A high degree of day-to-day variability in SARS-CoV-2 test collection was observed, consistent with uneven implementation of serial point prevalence surveys. Nonetheless, we were able to achieve excellent model calibration and estimate the basic reproductive number (R0) in all included facilities. The stochastic model captured differences between facilities; the estimated facility R0 ranged from 1.4 to 2.5 (Figure 2). We did not find evidence that a delay in collection of point prevalence survey was associated with an increased R0.ConclusionFitting simple dynamic models to COVID surveillance data enables a comparison of nursing homes with respect to facility R0 and an evaluation of the effect of different control policies. More complex models are needed to yield insights about the relative importance of alternative transmission pathways, such as combinations between staff and residents.DisclosuresAll Authors: No reported disclosures

Bayesian inference for nonlinear diffusions, observed at discrete times, is a challenging task that has prompted the development of a number of algorithms, mainly within the computational statistics community. We propose a new direction, and accompanying methodology—borrowing ideas from statistical physics and computational chemistry—for inferring the posterior distribution of latent diffusion paths and model parameters, given observations of the process. Joint configurations of the underlying process noise and of parameters, mapping onto diffusion paths consistent with observations, form an implicitly defined manifold. Then, by making use of a constrained Hamiltonian Monte Carlo algorithm on the embedded manifold, we are able to perform computationally efficient inference for a class of discretely observed diffusion models. Critically, in contrast with other approaches proposed in the literature, our methodology is highly automated, requiring minimal user intervention and applying alike in a range of settings, including: elliptic or hypo-elliptic systems; observations with or without noise; linear or non-linear observation operators. Exploiting Markovianity, we propose a variant of the method with complexity that scales linearly in the resolution of path discretisation and the number of observation times. Python code reproducing the results is available at http://doi.org/10.5281/zenodo.5796148.