Metropolis-Hastings Sampling Research Articles

A Heuristic Method for Power Pylon Reconstruction from Airborne LiDAR Data

Object reconstruction from airborne LiDAR data is a hot topic in photogrammetry and remote sensing. Power fundamental infrastructure monitoring plays a vital role in power transmission safety. This paper proposes a heuristic reconstruction method for power pylons widely used in high voltage transmission systems from airborne LiDAR point cloud, which combines both data-driven and model-driven strategies. Structurally, a power pylon can be decomposed into two parts: the pylon body and head. The reconstruction procedure assembles two parts sequentially: firstly, the pylon body is reconstructed by a data-driven strategy, where a RANSAC-based algorithm is adopted to fit four principal legs; secondly, a model-driven strategy is used to reconstruct the pylon head with the aid of a predefined 3D head model library, where the pylon head’s type is recognized by a shape context algorithm, and their parameters are estimated by a Metropolis–Hastings sampler coupled with a Simulated annealing algorithm. The proposed method has two advantages: (1) optimal strategies are adopted to reconstruct different pylon parts, which are robust to noise and partially missing data; and (2) both the number of parameters and their search space are greatly reduced when estimating the head model’s parameters, as the body reconstruction results information about the original point cloud, and relationships between parameters are used in the pylon head reconstruction process. Experimental results show that the proposed method can efficiently reconstruct power pylons, and the average residual between the reconstructed models and the raw data was smaller than 0.3 m.

Modelling and Computation Using NCoRM Mixtures for Density Regression

Normalized compound random measures are flexible nonparametric priors for related distributions. We consider building general nonparametric regression models using normalized compound random measure mixture models. Posterior inference is made using a novel pseudo-marginal Metropolis-Hastings sampler for normalized compound random measure mixture models. The algorithm makes use of a new general approach to the unbiased estimation of Laplace functionals of compound random measures (which includes completely random measures as a special case). The approach is illustrated on problems of density regression.

Nautilus

Natural images often exhibit symmetries that should be taken into account when editing them. In this paper we present Nautilus --- a method for automatically identifying symmetric regions in an image along with their corresponding symmetry transformations. We compute dense local similarity symmetry transformations using a novel variant of the Generalised PatchMatch algorithm that uses Metropolis-Hastings sampling. We combine and refine these local symmetries using an extended Lucas-Kanade algorithm to compute regional transformations and their spatial extents. Our approach produces dense estimates of complex symmetries that are combinations of translation, rotation, scale, and reflection under perspective distortion. This enables a number of automatic symmetry-aware image editing applications including inpainting, rectification, beautification, and segmentation, and we demonstrate state-of-the-art applications for each of them.

Pseudo-marginal Metropolis–Hastings sampling using averages of unbiased estimators

SummaryWe consider a pseudo-marginal Metropolis–Hastings kernel ${\mathbb{P}}_m$ that is constructed using an average of $m$ exchangeable random variables, and an analogous kernel ${\mathbb{P}}_s$ that averages $s<m$ of these same random variables. Using an embedding technique to facilitate comparisons, we provide a lower bound for the asymptotic variance of any ergodic average associated with ${\mathbb{P}}_m$ in terms of the asymptotic variance of the corresponding ergodic average associated with ${\mathbb{P}}_s$. We show that the bound is tight and disprove a conjecture that when the random variables to be averaged are independent, the asymptotic variance under ${\mathbb{P}}_m$ is never less than $s/m$ times the variance under ${\mathbb{P}}_s$. The conjecture does, however, hold for continuous-time Markov chains. These results imply that if the computational cost of the algorithm is proportional to $m$, it is often better to set $m=1$. We provide intuition as to why these findings differ so markedly from recent results for pseudo-marginal kernels employing particle filter approximations. Our results are exemplified through two simulation studies; in the first the computational cost is effectively proportional to $m$ and in the second there is a considerable start-up cost at each iteration.

Monopulse joint parameter estimation of multiple unresolved targets within the radar beam

Aiming at the problem of the parameter estimation of multiple unresolved targets within the radar beam, using the joint bin processing model, a method of jointly estimating the number and the position of the targets is proposed based on reversible jump Markov Chain Monte Carlo (RJ-MCMC). Reasonable assumptions of the prior distributions and Bayesian theory are adopted to obtain the posterior probability density function of the estimated parameters from the conditional likelihood function of the observation, and then the acceptance ratios of the birth, death and update moves are given. During the update move, a hybrid Metropolis-Hastings (MH) sampling algorithm is used to make a better exploration of the parameter space. The simulation results show that this new method outperforms the method of ML-MLD [11] proposed by X.Zhang for similar estimation accuracy is achieved while fewer sub-pulses are needed.

BayesWHAM: A Bayesian approach for free energy estimation, reweighting, and uncertainty quantification in the weighted histogram analysis method.

The weighted histogram analysis method (WHAM) is a powerful approach to estimate molecular free energy surfaces (FES) from biased simulation data. Bayesian reformulations of WHAM are valuable in proving statistically optimal use of the data and providing a transparent means to incorporate regularizing priors and estimate statistical uncertainties. In this work, we develop a fully Bayesian treatment of WHAM to generate statistically optimal FES estimates in any number of biasing dimensions under arbitrary choices of the Bayes prior. Rigorous uncertainty estimates are generated by Metropolis-Hastings sampling from the Bayes posterior. We also report a means to project the FES and its uncertainties into arbitrary auxiliary order parameters beyond those in which biased sampling was conducted. We demonstrate the approaches in applications of alanine dipeptide and the unthreading of a synthetic mimic of the astexin-3 lasso peptide. Open-source MATLAB and Python implementations of our codes are available for free public download. © 2017 Wiley Periodicals, Inc.

Model selection for profile structure in Bayesian geoacoustic inversion

This paper considers model selection in Bayesian geoacoustic inversion, specifically, the role of seabed parameterization in resolving geoacoustic profile structure. Bayesian inversion is formulated in terms of the posterior probability density (PPD) over the model parameters which are sampled numerically: Metropolis-Hastings sampling in principal-component space enhanced by parallel tempering is employed here. A key aspect of quantitative geoacoustic inversion is that of parameterizing the seabed model. Trans-dimensional (trans-D) inversion methods model the seabed as a sequence of discontinuous uniform layers and sample probabilistically over the number of layers. However, in some cases it may be expected for seabed properties to vary as smooth, continuous gradients which are not well represented by uniform layers. Most gradient-based inversions assume a representative functional form, such as a power law. A recent alternative is based on a linear combination of Bernstein-polynomial basis functions. This approach is more general and allows the form of the profile to be determined by the data, rather than by a subjective model choice. This paper compares trans-D, power-law, and Bernstein-polynomial inversions for the problem of estimating seabed shear-wave speed profiles from the dispersion of interface waves. Simulations and data from Oslofjorden and/or the North Sea will be considered.

Programming With Models: Writing Statistical Algorithms for General Model Structures With NIMBLE

ABSTRACTWe describe NIMBLE, a system for programming statistical algorithms for general model structures within R. NIMBLE is designed to meet three challenges: flexible model specification, a language for programming algorithms that can use different models, and a balance between high-level programmability and execution efficiency. For model specification, NIMBLE extends the BUGS language and creates model objects, which can manipulate variables, calculate log probability values, generate simulations, and query the relationships among variables. For algorithm programming, NIMBLE provides functions that operate with model objects using two stages of evaluation. The first stage allows specialization of a function to a particular model and/or nodes, such as creating a Metropolis-Hastings sampler for a particular block of nodes. The second stage allows repeated execution of computations using the results of the first stage. To achieve efficient second-stage computation, NIMBLE compiles models and functions via C++, using the Eigen library for linear algebra, and provides the user with an interface to compiled objects. The NIMBLE language represents a compilable domain-specific language (DSL) embedded within R. This article provides an overview of the design and rationale for NIMBLE along with illustrative examples including importance sampling, Markov chain Monte Carlo (MCMC) and Monte Carlo expectation maximization (MCEM). Supplementary materials for this article are available online.

Semistochastic Heat-Bath Configuration Interaction Method: Selected Configuration Interaction with Semistochastic Perturbation Theory.

We extend the recently proposed heat-bath configuration interaction (HCI) method [Holmes, Tubman, Umrigar, J. Chem. Theory Comput. 2016, 12, 3674], by introducing a semistochastic algorithm for performing multireference Epstein-Nesbet perturbation theory, in order to completely eliminate the severe memory bottleneck of the original method. The proposed algorithm has several attractive features. First, there is no sign problem that plagues several quantum Monte Carlo methods. Second, instead of using Metropolis-Hastings sampling, we use the Alias method to directly sample determinants from the reference wave function, thus avoiding correlations between consecutive samples. Third, in addition to removing the memory bottleneck, semistochastic HCI (SHCI) is faster than the deterministic variant for many systems if a stochastic error of 0.1 mHa is acceptable. Fourth, within the SHCI algorithm one can trade memory for a modest increase in computer time. Fifth, the perturbative calculation is embarrassingly parallel. The SHCI algorithm extends the range of applicability of the original algorithm, allowing us to calculate the correlation energy of very large active spaces. We demonstrate this by performing calculations on several first row dimers including F2 with an active space of (14e, 108o), Mn-Salen cluster with an active space of (28e, 22o), and Cr2 dimer with up to a quadruple-ζ basis set with an active space of (12e, 190o). For these systems we were able to obtain better than 1 mHa accuracy with a wall time of merely 55 s, 37 s, and 56 min on 1, 1, and 4 nodes, respectively.

Bayesian estimation of discretely observed multi-dimensional diffusion processes using guided proposals

Estimation of parameters of a diffusion based on discrete time observations poses a difficult problem due to the lack of a closed form expression for the likelihood. From a Bayesian computational perspective it can be casted as a missing data problem where the diffusion bridges in between discrete-time observations are missing. The computational problem can then be dealt with using a Markov-chain Monte-Carlo method known as data-augmentation. If unknown parameters appear in the diffusion coefficient, direct implementation of data-augmentation results in a Markov chain that is reducible. Furthermore, data-augmentation requires efficient sampling of diffusion bridges, which can be difficult, especially in the multidimensional case. We present a general framework to deal with with these problems that does not rely on discretisation. The construction generalises previous approaches and sheds light on the assumptions necessary to make these approaches work. We define a random-walk type Metropolis-Hastings sampler for updating diffusion bridges. Our methods are illustrated using guided proposals for sampling diffusion bridges. These are Markov processes obtained by adding a guiding term to the drift of the diffusion. We give general guidelines on the construction of these proposals and introduce a time change and scaling of the guided proposal that reduces discretisation error. Numerical examples demonstrate the performance of our methods.

Efficient Markov Chain Monte Carlo for combined Subset Simulation and nonlinear finite element analysis

Typical probabilistic problems in an engineering context include rare event probability estimation for physical models where spatial autocorrelation of material property parameters is significant. Subset Simulation, a Markov Chain Monte Carlo technique, can be used to estimate rare event probabilities in physical models more efficiently than Monte Carlo Simulation. This efficiency gain is important when the sampling operation is computationally demanding, as is the case in the solution of stochastic Partial Differential Equations. In high dimensional spaces where Polynomial Chaos or other direct integration techniques become intractable, sampling methods may be the only way to compute integral functions in probabilistic analysis. In this paper, Subset Simulation is applied to probability of failure estimation in nonlinear elasto-plastic finite element problems. Further, a derivation of confidence intervals for Subset Simulation relative errors is presented. This new technique allows for vastly improved efficiency in the computation of error estimates for Subset Simulation. Significantly, the numerical studies presented indicate that for the tested finite element problems, Metropolis–Hastings sampling can outperform Componentwise Metropolis–Hastings and Gibbs sampling. This result is relevant to the design of efficient Subset Simulation methodologies.

Not every Gibbs sampler is a special case of the Metropolis–Hastings algorithm

ABSTRACTIt is commonly asserted that the Gibbs sampler is a special case of the Metropolis–Hastings (MH) algorithm. While this statement is true for certain Gibbs samplers, it is not true in general for the version that is taught and used most often, namely, the deterministic scan Gibbs sampler. In this note, I prove that that there exist deterministic scan Gibbs samplers that do not exhibit detailed balance and hence cannot be considered MH samplers. The nuances of various Gibbs sampling schemes are discussed.

Uniform Sampling of Directed and Undirected Graphs Conditional on Vertex Connectivity

Many applications in graph analysis require a space of graphs or networks to be sampled uniformly at random. For example, one may need to efficiently draw a small representative sample of graphs from a particular large target space. We assume that a uniform distribution f(N,E)=1/|X| has been defined, where N is a set of nodes, E is a set of edges, (N, E) is a graph in the target space X and |X| is the (finite) total number of graphs in the target space. We propose a new approach to sample graphs at random from such a distribution. The new approach uses a Markov chain Monte Carlo method called the Neighbourhood Sampler. We validate the new sampling technique by simulating from feasible spaces of directed or undirected graphs, and compare its computational efficiency with the conventional Metropolis-Hastings Sampler. The simulation results indicate efficient uniform sampling of the target spaces, and more rapid rate of convergence than Metropolis-Hastings Sampler.

Cure Models based on Weibull Distribution with and without Covariates using Right Censored Data

In this paper we use a methodology based on the Weibull distributions covariates in the presence of cure fraction models, censored data and covariates. Objective: The objective of the study is to check the performance of mixture and non-mixture cure models based on LPML. Methods/Analysis: Two models were explored here in which are the mixture and non-mixture cure fraction models. Inferences for the models are obtained under the Bayesian approach via Markov Chain Monte Carlo (MCMC) where the posterior estimates were obtained by using Metropolis-Hastings sampling methods in the presence of covariates and without covariates considering a real life time dataset and comparing the two cure models using the Log Pseudo Maximum Likelihood estimates (LPML) and some related special cases of the distribution. Findings/ Conclusion: We observed that the Weibull distribution has the least LPML value while its special cases where the two models are quite similar having the highest values on the other hand, the Mixture fits better than the non-mixture having the highest (LPML) based on the results obtain from all the models suggesting that the standard parametric cure (mixture) model fits the AML data which shows a great indication of similarity with the covariates and flexibility of the models.

Granger Causality and Regime Inference in Markov Switching VAR Models with Bayesian Methods

Summary In this paper, we derive restrictions for Granger noncausality in MS-VAR models and show under what conditions a variable does not affect the forecast of the hidden Markov process. To assess the noncausality hypotheses, we apply Bayesian inference. The computational tools include a novel block Metropolis–Hastings sampling algorithm for the estimation of the underlying models. We analyze a system of monthly US data on money and income. The results of testing in MS-VARs contradict those obtained with linear VARs: the money aggregate M1 helps in forecasting industrial production and in predicting the next period's state. Copyright © 2016 John Wiley & Sons, Ltd.

WarpLDA

Developing efficient and scalable algorithms for Latent Dirichlet Allocation (LDA) is of wide interest for many applications. Previous work has developed an O (1) Metropolis-Hastings (MH) sampling method for each token. However, its performance is far from being optimal due to frequent cache misses caused by random accesses to the parameter matrices. In this paper, we first carefully analyze the memory access behavior of existing algorithms for LDA by cache locality at document level. We then develop WarpLDA, which achieves O (1) time complexity per-token and fits the randomly accessed memory perdocument in the L3 cache. Our empirical results in a wide range of testing conditions demonstrate that WarpLDA is consistently 5-15x faster than the state-of-the-art MH-based LightLDA, and is faster than the state-of-the-art sparsity aware F+LDA in most settings. Our WarpLDA learns a million topics from 639 millions of documents in only five hours at an unprecedented throughput of 11 billion tokens per second.

A BAYESIAN HIERARCHICAL SPATIAL MODEL FOR DENTAL CARIES ASSESSMENT USING NON-GAUSSIAN MARKOV RANDOM FIELDS.

Research in dental caries generates data with two levels of hierarchy: that of a tooth overall and that of the different surfaces of the tooth. The outcomes often exhibit spatial referencing among neighboring teeth and surfaces, i.e., the disease status of a tooth or surface might be influenced by the status of a set of proximal teeth/surfaces. Assessments of dental caries (tooth decay) at the tooth level yield binary outcomes indicating the presence/absence of teeth, and trinary outcomes at the surface level indicating healthy, decayed, or filled surfaces. The presence of these mixed discrete responses complicates the data analysis under a unified framework. To mitigate complications, we develop a Bayesian two-level hierarchical model under suitable (spatial) Markov random field assumptions that accommodates the natural hierarchy within the mixed responses. At the first level, we utilize an autologistic model to accommodate the spatial dependence for the tooth-level binary outcomes. For the second level and conditioned on a tooth being non-missing, we utilize a Potts model to accommodate the spatial referencing for the surface-level trinary outcomes. The regression models at both levels were controlled for plausible covariates (risk factors) of caries, and remain connected through shared parameters. To tackle the computational challenges in our Bayesian estimation scheme caused due to the doubly-intractable normalizing constant, we employ a double Metropolis-Hastings sampler. We compare and contrast our model performances to the standard non-spatial (naive) model using a small simulation study, and illustrate via an application to a clinical dataset on dental caries.

Hamiltonian Monte Carlo Sampling in Bayesian Empirical Likelihood Computation

SummaryWe consider Bayesian empirical likelihood estimation and develop an efficient Hamiltonian Monte Carlo method for sampling from the posterior distribution of the parameters of interest. The method proposed uses hitherto unknown properties of the gradient of the underlying log-empirical-likelihood function. We use results from convex analysis to show that these properties hold under minimal assumptions on the parameter space, prior density and the functions used in the estimating equations determining the empirical likelihood. Our method employs a finite number of estimating equations and observations but produces valid semiparametric inference for a large class of statistical models including mixed effects models, generalized linear models and hierarchical Bayes models. We overcome major challenges posed by complex, non-convex boundaries of the support routinely observed for empirical likelihood which prevent efficient implementation of traditional Markov chain Monte Carlo methods like random-walk Metropolis–Hastings sampling etc. with or without parallel tempering. A simulation study confirms that our method converges quickly and draws samples from the posterior support efficiently. We further illustrate its utility through an analysis of a discrete data set in small area estimation.

Inferences on target speed, depth, and range from a continuous wave transmission

Tracking an underwater mobile object by means of a continuous wave transmission is challenging in part due to the difficulty of drawing inferences on closely spaced tones associated with target depth and speed. Considered here is a bistatic sonar arrangement employed to infer the depth, speed and range of an oncoming submerged object. Computation of the full posterior probability distribution of the returned amplitudes and frequencies from both the prior distribution and a finite duration window of the received waveform is made by Markov-chain Monte Carlo sampling. A Gibbs sampler is employed to construct the posterior joint density of all parameters by taking full advantage of the analytic tractability of the conditional and marginal densities of the received amplitudes while those of the ordered frequencies are constructed numerically by either an inverse quantile sampling or a Metropolis-Hastings sampler. The inferred density of depth, range, and speed of the target is accomplished by constructing a numerical inverse-transformation of the forward propagation model.

Bayesian localization, tracking, and environmental inference

This paper presents an overview of a general Bayesian inference approach to source localization, tracking and/or environmental estimation. Source location and spectral parameters together with environmental parameters are all considered unknown random variables to be estimated from prior information and observed acoustic data. The relative level of prior information for various parameters differentiates applications of interest. For instance, controlled-source geoacoustic inversion typically involves large prior uncertainties for seabed parameters but small uncertainties for source locations, although some applications, such as inverting noise from ships-of-opportunity, may involve larger location uncertainties. Alternatively, source localization in an uncertain environment typically involves non-informative location priors and environmental priors that reflect available knowledge. Tracking applications include additional prior constraints on source speed. In all cases, the goal is to compute marginals of the posterior probability density for source and environmental parameters, quantifying the information content of the data and prior. This is typically carried out with Markov-chain Monte Carlo methods including Metropolis-Hastings sampling and/or Gibbs sampling, with various approaches applied to improve efficiency (e.g., principal-component sampling, parallel tempering) and generality (trans-dimensional inversion). Multiple-source localization minimizes the Bayesian information criterion to estimate the number of sources.