Dynamically Rescaled Hamiltonian Monte Carlo for Bayesian Hierarchical Models

Dynamically rescaled Hamiltonian Monte Carlo is introduced as a computationally fast and easily implemented method for performing full Bayesian analysis in hierarchical statistical models. The method relies on introducing a modified parameterization so that the reparameterized target distribution has close to constant scaling properties, and thus is easily sampled using standard (Euclidian metric) Hamiltonian Monte Carlo. Provided that the parameterizations of the conditional distributions specifying the hierarchical model are “constant information parameterizations” (CIPs), the relation between the modified- and original parameterization is bijective, explicitly computed, and admit exploitation of sparsity in the numerical linear algebra involved. CIPs for a large catalogue of statistical models are presented, and from the catalogue, it is clear that many CIPs are currently routinely used in statistical computing. A relation between the proposed methodology and a class of explicitly integrated Riemann manifold Hamiltonian Monte Carlo methods is discussed. The methodology is illustrated on several example models, including a model for inflation rates with multiple levels of nonlinearly dependent latent variables. Supplementary materials for this article are available online.     

Noisy Hamiltonian Monte Carlo for Doubly Intractable Distributions

Hamiltonian Monte Carlo (HMC) has been progressively incorporated within the statistician’s toolbox as an alternative sampling method in settings when standard Metropolis–Hastings is inefficient. HMC generates a Markov chain on an augmented state space with transitions based on a deterministic differential flow derived from Hamiltonian mechanics. In practice, the evolution of Hamiltonian systems cannot be solved analytically, requiring numerical integration schemes. Under numerical integration, the resulting approximate solution no longer preserves the measure of the target distribution, therefore an accept–reject step is used to correct the bias. For doubly intractable distributions—such as posterior distributions based on Gibbs random fields—HMC suffers from some computational difficulties: computation of gradients in the differential flow and computation of the accept–reject proposals poses difficulty. In this article, we study the behavior of HMC when these quantities are replaced by Monte Carlo estimates. Supplemental codes for implementing methods used in the article are available online.     

基于贝叶斯MCMC算法的美式期权定价

鉴于美式期权的定价具有后向迭代搜索特征,本文结合Longstaff和Schwartz提出的美式期权定价的最小二乘模拟方法,研究基于马尔科夫链蒙特卡洛算法对回归方程系数的估计,实现对美式期权的双重模拟定价.通过对无红利美式看跌股票期权定价进行大量实证模拟,从期权价值定价误差等方面同著名的最小二乘蒙特卡洛模拟方法进行对比分析,结果表明基于MCMC回归算法给出的美式期权定价具有更高的精确度.模拟实证结果表明本文提出的对美式期权定价方法具有较好的可行性、有效性与广泛的适用性.该方法的不足之处就是类似于一般的蒙特卡洛方法,会使得求解的计算量有所加大.     

Markov Chain Monte Carlo from Lagrangian Dynamics

Hamiltonian Monte Carlo (HMC) improves the computational efficiency of the Metropolis–Hastings algorithm by reducing its random walk behavior. Riemannian HMC (RHMC) further improves the performance of HMC by exploiting the geometric properties of the parameter space. However, the geometric integrator used for RHMC involves implicit equations that require fixed-point iterations. In some cases, the computational overhead for solving implicit equations undermines RHMC’s benefits. In an attempt to circumvent this problem, we propose an explicit integrator that replaces the momentum variable in RHMC by velocity. We show that the resulting transformation is equivalent to transforming Riemannian Hamiltonian dynamics to Lagrangian dynamics. Experimental results suggest that our method improves RHMC’s overall computational efficiency in the cases considered. All computer programs and datasets are available online (http://www.ics.uci.edu/babaks/Site/Codes.html) to allow replication of the results reported in this article.     

Hybrid Monte Carlo on Hilbert spaces

The Hybrid Monte Carlo (HMC) algorithm provides a framework for sampling from complex, high-dimensional target distributions. In contrast with standard Markov chain Monte Carlo (MCMC) algorithms, it generates nonlocal, nonsymmetric moves in the state space, alleviating random walk type behaviour for the simulated trajectories. However, similarly to algorithms based on random walk or Langevin proposals, the number of steps required to explore the target distribution typically grows with the dimension of the state space. We define a generalized HMC algorithm which overcomes this problem for target measures arising as finite-dimensional approximations of measures π which have density with respect to a Gaussian measure on an infinite-dimensional Hilbert space. The key idea is to construct an MCMC method which is well defined on the Hilbert space itself.We successively address the following issues in the infinite-dimensional setting of a Hilbert space: (i) construction of a probability measure Π in an enlarged phase space having the target π as a marginal, together with a Hamiltonian flow that preserves Π; (ii) development of a suitable geometric numerical integrator for the Hamiltonian flow; and (iii) derivation of an accept/reject rule to ensure preservation of Π when using the above numerical integrator instead of the actual Hamiltonian flow. Experiments are reported that compare the new algorithm with standard HMC and with a version of the Langevin MCMC method defined on a Hilbert space.     

蒙特卡罗模拟法在边坡可靠性分析中的运用

江永红．蒙特卡罗模拟法在边坡可靠性分析中的运用．数理统计与管理，１９９８，１７（１），１３～１６．可靠性分析是边坡工程及滑坡治理中的重要研究课题。鉴于决定边坡可靠性的诸变量多为随机变量，本文论述了用蒙特卡罗模拟法计算边坡可靠度的基本原理，对模拟次数确定、计算误差估计等问题提出了解决办法，并结合具体运用说明该方法的实施步骤     

On Markov chain Monte Carlo Algorithms for Computing Conditional Expectations Based on Sufficient Statistics

Much work has focused on developing exact tests for the analysis of discrete data using log linear or logistic regression models. A parametric model is tested for a dataset by conditioning on the value of a sufficient statistic and determining the probability of obtaining another dataset as extreme or more extreme relative to the general model, where extremeness is determined by the value of a test statistic such as the chi-square or the log-likelihood ratio. Exact determination of these probabilities can be infeasible for high dimensional problems, and asymptotic approximations to them are often inaccurate when there are small data entries and/or there are many nuisance parameters. In these cases Monte Carlo methods can be used to estimate exact probabilities by randomly generating datasets (tables) that match the sufficient statistic of the original table. However, naive Monte Carlo methods produce tables that are usually far from matching the sufficient statistic. The Markov chain Monte Carlo method used in this work (the regression/attraction approach) uses attraction to concentrate the distribution around the set of tables that match the sufficient statistic, and uses regression to take advantage of information in tables that “almost” match. It is also more general than others in that it does not require the sufficient statistic to be linear, and it can be adapted to problems involving continuous variables. The method is applied to several high dimensional settings including four-way tables with a model of no four-way interaction, and a table of continuous data based on beta distributions. It is powerful enough to deal with the difficult problem of four-way tables and flexible enough to handle continuous data with a nonlinear sufficient statistic.     

Monte Carlo optimization

Monte Carlo optimization techniques for solving mathematical programming problems have been the focus of some debate. This note reviews the debate and puts these stochastic methods in their proper perspective.     

The Monte Carlo EM method for estimating multinomial probit latent variable models

We propose a multinomial probit (MNP) model that is defined by a factor analysis model with covariates for analyzing unordered categorical data, and discuss its identification. Some useful MNP models are special cases of the proposed model. To obtain maximum likelihood estimates, we use the EM algorithm with its M-step greatly simplified under Conditional Maximization and its E-step made feasible by Monte Carlo simulation. Standard errors are calculated by inverting a Monte Carlo approximation of the information matrix using Louis’s method. The methodology is illustrated with a simulated data.     

Computing the variance of a conditional expectation via non-nested Monte Carlo

Computing the variance of a conditional expectation has often been of importance in uncertainty quantification. Sun et al. has introduced an unbiased nested Monte Carlo estimator, which they call $1\frac{1}{2}$-level simulation since the optimal inner-level sample size is bounded as the computational budget increases. In this letter, we construct unbiased non-nested Monte Carlo estimators based on the so-called pick-freeze scheme due to Sobol’. An extension of our approach to compute higher order moments of a conditional expectation is also discussed.     

Unbiased multi-index Monte Carlo

We introduce a new class of Monte Carlo-based approximations of expectations of random variables such that their laws are only available via certain discretizations. Sampling from the discretized versions of these laws can typically introduce a bias. In this paper, we show how to remove that bias, by introducing a new version of multi-index Monte Carlo (MIMC) that has the added advantage of reducing the computational effort, relative to i.i.d. sampling from the most precise discretization, for a given level of error. We cover extensions of results regarding variance and optimality criteria for the new approach. We apply the methodology to the problem of computing an unbiased mollified version of the solution of a partial differential equation with random coefficients. A second application concerns the Bayesian inference (the smoothing problem) of an infinite-dimensional signal modeled by the solution of a stochastic partial differential equation that is observed on a discrete space grid and at discrete times. Both applications are complemented by numerical simulations.     

Monte Carlo EM加速算法

EM算法是近年来常用的求后验众数的估计的一种数据增广算法, 但由于求出其E步中积分的显示表达式有时很困难, 甚至不可能, 限制了其应用的广泛性. 而Monte Carlo EM算法很好地解决了这个问题, 将EM算法中E步的积分用Monte Carlo模拟来有效实现, 使其适用性大大增强. 但无论是EM算法, 还是Monte Carlo EM算法, 其收敛速度都是线性的, 被缺损信息的倒数所控制, 当缺损数据的比例很高时, 收敛速度就非常缓慢. 而Newton-Raphson算法在后验众数的附近具有二次收敛速率. 本文提出Monte Carlo EM加速算法, 将Monte Carlo EM算法与Newton-Raphson算法结合, 既使得EM算法中的E步用Monte Carlo模拟得以实现, 又证明了该算法在后验众数附近具有二次收敛速度. 从而使其保留了Monte Carlo EM算法的优点, 并改进了Monte Carlo EM算法的收敛速度. 本文通过数值例子, 将Monte Carlo EM加速算法的结果与EM算法、Monte Carlo EM算法的结果进行比较, 进一步说明了Monte Carlo EM加速算法的优良性.     

多级运载火箭动态级间分离的Monte Carlo仿真

提出了研究多级运载火箭级间冷、热分离的计算公式．利用Monte Carlo仿真,考虑箭体分离时,采用远离正常设计值的参数去评估分离失败的风险．分析了各种扰动,如动态不平衡,剩余推力,分离机构引起的分离扰动,以及冷、热分离未对准的影响,得到分离箭体不会发生碰撞的条件．结果表明,当前设计满足分离要求．     

Useful Monte Carlo optimization

This contribution to the debate on Monte Carlo optimization methods shows that there exist techniques that may be useful in many technical applications.     

Multilevel sequential Monte Carlo: Mean square error bounds under verifiable conditions

In this article, we consider the multilevel sequential Monte Carlo (MLSMC) method of Beskos et al. (Stoch. Proc. Appl. [to appear]). This is a technique designed to approximate expectations w.r.t. probability laws associated to a discretization. For instance, in the context of inverse problems, where one discretizes the solution of a partial differential equation. The MLSMC approach is especially useful when independent, coupled sampling is not possible. Beskos et al. show that for MLSMC the computational effort to achieve a given error, can be less than independent sampling. In this article we significantly weaken the assumptions of Beskos et al., extending the proofs to non-compact state-spaces. The assumptions are based upon multiplicative drift conditions as in Kontoyiannis and Meyn (Electron. J. Probab. 10 [2005]: 61–123). The assumptions are verified for an example.     

Conditioning in Markov Chain Monte Carlo

Abstract

The so-called “Rao-Blackwellized” estimators proposed by Gelfand and Smith do not always reduce variance in Markov chain Monte Carlo when the dependence in the Markov chain is taken into account. An illustrative example is given, and a theorem characterizing the necessary and sufficient condition for such an estimator to always reduce variance is proved.     

Combinatorics of reliability Monte Carlo

We study one class of unbiased Monte Carlo estimators for system reliability, avoiding the rare event difficulty. This class is closely related to the system combinatorics and contains unique “extreme” members, having the minimum and maximum relative error. Some known Monte Carlo heuristics for network reliability, including fully polynomial cases, are of this type. ©1999 John Wiley & Sons, Inc. Random Struct. Alg., 14: 329–343, 1999     

Sequential Monte Carlo Methods for Option Pricing

In this article, we provide a review and development of sequential Monte Carlo (SMC) methods for option pricing. SMC are a class of Monte Carlo-based algorithms, that are designed to approximate expectations w.r.t a sequence of related probability measures. These approaches have been used successfully for a wide class of applications in engineering, statistics, physics, and operations research. SMC methods are highly suited to many option pricing problems and sensitivity/Greek calculations due to the nature of the sequential simulation. However, it is seldom the case that such ideas are explicitly used in the option pricing literature. This article provides an up-to-date review of SMC methods, which are appropriate for option pricing. In addition, it is illustrated how a number of existing approaches for option pricing can be enhanced via SMC. Specifically, when pricing the arithmetic Asian option w.r.t a complex stochastic volatility model, it is shown that SMC methods provide additional strategies to improve estimation.     

Bayesian estimation of NIG models via Markov chain Monte Carlo methods

The normal inverse Gaussian (NIG) distribution is a promising alternative for modelling financial data since it is a continuous distribution that allows for skewness and fat tails. There is an increasing number of applications of the NIG distribution to financial problems. Due to the complicated nature of its density, estimation procedures are not simple. In this paper we propose Bayesian estimation for the parameters of the NIG distribution via an MCMC scheme based on the Gibbs sampler. Our approach makes use of the data augmentation provided by the mixture representation of the distribution. We also extend the model to allow for modelling heteroscedastic regression situations. Examples with financial and simulated data are provided. Copyright © 2004 John Wiley & Sons, Ltd.     

Partition-Weighted Monte Carlo Estimation

Although various efficient and sophisticated Markov chain Monte Carlo sampling methods have been developed during the last decade, the sample mean is still a dominant in computing Bayesian posterior quantities. The sample mean is simple, but may not be efficient. The weighted sample mean is a natural generalization of the sample mean. In this paper, a new weighted sample mean is proposed by partitioning the support of posterior distribution, so that the same weight is assigned to observations that belong to the same subset in the partition. A novel application of this new weighted sample mean in computing ratios of normalizing constants and necessary theory are provided. Illustrative examples are given to demonstrate the methodology.