Optimal Scaling of Random Walk Metropolis Algorithms with Non-Gaussian Proposals

The asymptotic optimal scaling of random walk Metropolis (RWM) algorithms with Gaussian proposal distributions is well understood for certain specific classes of target distributions. These asymptotic results easily extend to any light tailed proposal distribution with finite fourth moment. However, heavy tailed proposal distributions such as the Cauchy distribution are known to have a number of desirable properties, and in many situations are preferable to light tailed proposal distributions. Therefore we consider the question of scaling for Cauchy distributed proposals for a wide range of independent and identically distributed (iid) product densities. The results are somewhat surprising as to when and when not Cauchy distributed proposals are preferable to Gaussian proposal distributions. This provides motivation for finding proposal distributions which improve on both Gaussian and Cauchy proposals, both for finite dimensional target distributions and asymptotically as the dimension of the target density, d → ∞. Throughout we seek the scaling of the proposal distribution which maximizes the expected squared jumping distance (ESJD).     

Posterior Simulation in the Generalized Linear Mixed Model With Semiparametric Random Effects

Generalized linear mixed models with semiparametric random effects are useful in a wide variety of Bayesian applications. When the random effects arise from a mixture of Dirichlet process (MDP) model with normal base measure, Gibbs samplingalgorithms based on the Pólya urn scheme are often used to simulate posterior draws in conjugate models (essentially, linear regression models and models for binary outcomes). In the nonconjugate case, some common problems associated with existing simulation algorithms include convergence and mixing difficulties.

This article proposes an algorithm for MDP models with exponential family likelihoods and normal base measures. The algorithm proceeds by making a Laplace approximation to the likelihood function, thereby matching the proposal with that of the Gibbs sampler. The proposal is accepted or rejected via a Metropolis-Hastings step. For conjugate MDP models, the algorithm is identical to the Gibbs sampler. The performance of the technique is investigated using a Poisson regression model with semi-parametric random effects. The algorithm performs efficiently and reliably, even in problems where large-sample results do not guarantee the success of the Laplace approximation. This is demonstrated by a simulation study where most of the count data consist of small numbers. The technique is associated with substantial benefits relative to existing methods, both in terms of convergence properties and computational cost.     

Optimal Scaling for Random Walk Metropolis on Spherically Constrained Target Densities

We consider the problem of optimal scaling of the proposal variance for multidimensional random walk Metropolis algorithms. It is well known, for a wide range of continuous target densities, that the optimal scaling of the proposal variance leads to an average acceptance rate of 0.234. Therefore a natural question is, do similar results hold for target densities which have discontinuities? In the current work, we answer in the affirmative for a class of spherically constrained target densities. Even though the acceptance probability is more complicated than for continuous target densities, the optimal scaling of the proposal variance again leads to an average acceptance rate of 0.234.     

Efficient Adaptive MCMC Through Precision Estimation

The performance of Markov chain Monte Carlo (MCMC) algorithms like the Metropolis Hastings Random Walk (MHRW) is highly dependent on the choice of scaling matrix for the proposal distributions. A popular choice of scaling matrix in adaptive MCMC methods is to use the empirical covariance matrix (ECM) of previous samples. However, this choice is problematic if the dimension of the target distribution is large, since the ECM then converges slowly and is computationally expensive to use. We propose two algorithms to improve convergence and decrease computational cost of adaptive MCMC methods in cases when the precision (inverse covariance) matrix of the target density can be well-approximated by a sparse matrix. The first is an algorithm for online estimation of the Cholesky factor of a sparse precision matrix. The second estimates the sparsity structure of the precision matrix. Combining the two algorithms allows us to construct precision-based adaptive MCMC algorithms that can be used as black-box methods for densities with unknown dependency structures. We construct precision-based versions of the adaptive MHRW and the adaptive Metropolis adjusted Langevin algorithm and demonstrate the performance of the methods in two examples. Supplementary materials for this article are available online.     

Scaling Analysis of Delayed Rejection MCMC Methods

In this paper, we study the asymptotic efficiency of the delayed rejection strategy. In particular, the efficiency of the delayed rejection Metropolis–Hastings algorithm is compared to that of the regular Metropolis algorithm. To allow for a fair comparison, the study is carried under optimal mixing conditions for each of these algorithms. After introducing optimal scaling results for the delayed rejection (DR) algorithm, we outline the fact that the second proposal after the first rejection is discarded, with a probability tending to 1 as the dimension of the target density increases. To overcome this drawback, a modification of the delayed rejection algorithm is proposed, in which the direction of the different proposals is fixed once for all, and the Metropolis–Hastings accept-reject mechanism is used to select a proper scaling along the search direction. It is shown that this strategy significantly outperforms the original DR and Metropolis algorithms, especially when the dimension becomes large. We include numerical studies to validate these conclusions.     

Variable target value subgradient method

Polyak's subgradient algorithm for nondifferentiable optimization problems requires prior knowledge of the optimal value of the objective function to find an optimal solution. In this paper we extend the convergence properties of the Polyak's subgradient algorithm with a fixed target value to a more general case with variable target values. Then a target value updating scheme is provided which finds an optimal solution without prior knowledge of the optimal objective value. The convergence proof of the scheme is provided and computational results of the scheme are reported.Most of this research was performed when the first author was visiting Decision and Information Systems Department, College of Business, Arizona State University.     

Identification of Regeneration Times in MCMC Simulation,With Application to Adaptive Schemes

Regeneration is a useful tool in Markov chain Monte Carlo simulation because it can be used to side-step the burn-in problem and to construct better estimates of the variance of parameter estimates themselves. It also provides a simple way to introduce adaptive behavior into a Markov chain, and to use parallel processors to build a single chain. Regeneration is often difficult to take advantage of because, for most chains, no recurrent proper atom exists, and it is not always easy to use Nummelin's splitting method to identify regeneration times. This article describes a constructive method for generating a Markov chain with a specified target distribution and identifying regeneration times. As a special case of the method, an algorithm which can be “wrapped” around an existing Markov transition kernel is given. In addition, a specific rule for adapting the transition kernel at regeneration times is introduced, which gradually replaces the original transition kernel with an independence-sampling Metropolis-Hastings kernel using a mixture normal approximation to the target density as its proposal density. Computational gains for the regenerative adaptive algorithm are demonstrated in examples.     

Optimal Scaling for the Pseudo-Marginal Random Walk Metropolis: Insensitivity to the Noise Generating Mechanism

We examine the optimal scaling and the efficiency of the pseudo-marginal random walk Metropolis algorithm using a recently-derived result on the limiting efficiency as the dimension, \(d\rightarrow \infty \). We prove that the optimal scaling for a given target varies by less than 20 % across a wide range of distributions for the noise in the estimate of the target, and that any scaling that is within 20 % of the optimal one will be at least 70 % efficient. We demonstrate that this phenomenon occurs even outside the range of noise distributions for which we rigorously prove it. We then conduct a simulation study on an example with d = 10 where importance sampling is used to estimate the target density; we also examine results available from an existing simulation study with d = 5 and where a particle filter was used. Our key conclusions are found to hold in these examples also.     

On a Generalization of the Preconditioned Crank–Nicolson Metropolis Algorithm

Metropolis algorithms for approximate sampling of probability measures on infinite dimensional Hilbert spaces are considered, and a generalization of the preconditioned Crank–Nicolson (pCN) proposal is introduced. The new proposal is able to incorporate information on the measure of interest. A numerical simulation of a Bayesian inverse problem indicates that a Metropolis algorithm with such a proposal performs independently of the state-space dimension and the variance of the observational noise. Moreover, a qualitative convergence result is provided by a comparison argument for spectral gaps. In particular, it is shown that the generalization inherits geometric convergence from the Metropolis algorithm with pCN proposal.     

Implementing random scan Gibbs samplers

Summary  The Gibbs sampler, being a popular routine amongst Markov chain Monte Carlo sampling methodologies, has revolutionized the application of Monte Carlo methods in statistical computing practice. The performance of the Gibbs sampler relies heavily on the choice of sweep strategy, that is, the means by which the components or blocks of the random vector X of interest are visited and updated. We develop an automated, adaptive algorithm for implementing the optimal sweep strategy as the Gibbs sampler traverses the sample space. The decision rules through which this strategy is chosen are based on convergence properties of the induced chain and precision of statistical inferences drawn from the generated Monte Carlo samples. As part of the development, we analytically derive closed form expressions for the decision criteria of interest and present computationally feasible implementations of the adaptive random scan Gibbs sampler via a Gaussian approximation to the target distribution. We illustrate the results and algorithms presented by using the adaptive random scan Gibbs sampler developed to sample multivariate Gaussian target distributions, and screening test and image data. Research by RL and ZY supported in part by a US National Science Foundation FRG grant 0139948 and a grant from Lawrence Livermore National Laboratory, Livermore, California, USA.     

On the stability and ergodicity of adaptive scaling Metropolis algorithms

The stability and ergodicity properties of two adaptive random walk Metropolis algorithms are considered. Both algorithms adjust the scaling of the proposal distribution continuously based on the observed acceptance probability. Unlike the previously proposed forms of the algorithms, the adapted scaling parameter is not constrained within a predefined compact interval. The first algorithm is based on scale adaptation only, while the second one also incorporates covariance adaptation. A strong law of large numbers is shown to hold assuming that the target density is smooth enough and has either compact support or super-exponentially decaying tails.     

Gibbs Sampling Will Fail in Outlier Problems with Strong Masking

Abstract

This article discusses the convergence of the Gibbs sampling algorithm when it is applied to the problem of outlier detection in regression models. Given any vector of initial conditions, theoretically, the algorithm converges to the true posterior distribution. However, the speed of convergence may slow down in a high-dimensional parameter space where the parameters are highly correlated. We show that the effect of the leverage in regression models makes very difficult the convergence of the Gibbs sampling algorithm in sets of data with strong masking. The problem is illustrated with examples.     

Superlinearly convergent affine scaling interior trust-region method for linear constrained LC<Superscript>1</Superscript> minimization

We extend the classical affine scaling interior trust region algorithm for the linear constrained smooth minimization problem to the nonsmooth case where the gradient of objective function is only locally Lipschitzian. We propose and analyze a new affine scaling trust-region method in association with nonmonotonic interior backtracking line search technique for solving the linear constrained LC1 optimization where the second-order derivative of the objective function is explicitly required to be locally Lipschitzian. The general trust region subproblem in the proposed algorithm is defined by minimizing an augmented affine scaling quadratic model which requires both first and second order information of the objective function subject only to an affine scaling ellipsoidal constraint in a null subspace of the augmented equality constraints. The global convergence and fast local convergence rate of the proposed algorithm are established under some reasonable conditions where twice smoothness of the objective function is not required. Applications of the algorithm to some nonsmooth optimization problems are discussed.     

Stability of Merton's portfolio optimization problem for Lévy models

Merton's classical portfolio optimization problem for an investor, who can trade in a risk-free bond and a stock, can be extended to the case where the driving noise of the logreturns is a pure jump process instead of a Brownian motion. Benth et al. [4,5] solved the problem and found the optimal control implicitly given by an integral equation in the hyperbolic absolute risk aversion (HARA) utility case. There are several ways to approximate a Levy process with infinite activity by neglecting the small jumps or approximating them with a Brownian motion, as discussed in Asmussen and Rosinski [1]. In this setting, we study stability of the corresponding optimal investment problems. The optimal controls are solutions of integral equations, for which we study convergence. We are able to characterize the rate of convergence in terms of the variance of the small jumps. Additionally, we prove convergence of the corresponding wealth processes and indirect utilities (value functions).     

Real-Time Density and Mode Estimation With Application to Time-Dynamic Mode Tracking

We introduce a nonparametric time-dynamic kernel type density estimate for the situation where an underlying multivariate distribution evolves with time. Based on this time-dynamic density estimate, we propose nonparametric estimates for the time-dynamic mode of the underlying distribution. Our estimators involve boundary kernels for the time dimension so that the estimator is always centered at current time, and multivariate kernels for the spatial dimension of the time-evolving distribution. Under certain mild conditions, the asymptotic behavior of density and mode estimators, especially their uniform convergence in both time and space, is derived. A time-dynamic algorithm for mode tracking is proposed, including automatic bandwidth choices, and is implemented via a mean update algorithm. Simulation studies and real data illustrations demonstrate that the proposed methods work well in practice.     

Approximate Predetermined Convergence Properties of the Gibbs Sampler

This article aims to provide a method for approximately predetermining convergence properties of the Gibbs sampler. This is to be done by first finding an approximate rate of convergence for a normal approximation of the target distribution. The rates of convergence for different implementation strategies of the Gibbs sampler are compared to find the best one. In general, the limiting convergence properties of the Gibbs sampler on a sequence of target distributions (approaching a limit) are not the same as the convergence properties of the Gibbs sampler on the limiting target distribution. Theoretical results are given in this article to justify that under conditions, the convergence properties of the Gibbs sampler can be approximated as well. A number of practical examples are given for illustration.     

Probabilistic analysis of a differential equation for linear programming

In this paper we address the complexity of solving linear programming problems with a set of differential equations that converge to a fixed point that represents the optimal solution. Assuming a probabilistic model, where the inputs are i.i.d. Gaussian variables, we compute the distribution of the convergence rate to the attracting fixed point. Using the framework of Random Matrix Theory, we derive a simple expression for this distribution in the asymptotic limit of large problem size. In this limit, we find the surprising result that the distribution of the convergence rate is a scaling function of a single variable. This scaling variable combines the convergence rate with the problem size (i.e., the number of variables and the number of constraints). We also estimate numerically the distribution of the computation time to an approximate solution, which is the time required to reach a vicinity of the attracting fixed point. We find that it is also a scaling function. Using the problem size dependence of the distribution functions, we derive high probability bounds on the convergence rates and on the computation times to the approximate solution.     

First Hitting Time Analysis of the Independence Metropolis Sampler

In this paper, we study a special case of the Metropolis algorithm, the Independence Metropolis Sampler (IMS), in the finite state space case. The IMS is often used in designing components of more complex Markov Chain Monte Carlo algorithms. We present new results related to the first hitting time of individual states for the IMS. These results are expressed mostly in terms of the eigenvalues of the transition kernel. We derive a simple form formula for the mean first hitting time and we show tight lower and upper bounds on the mean first hitting time with the upper bound being the product of two factors: a “local” factor corresponding to the target state and a “global” factor, common to all the states, which is expressed in terms of the total variation distance between the target and the proposal probabilities. We also briefly discuss properties of the distribution of the first hitting time for the IMS and analyze its variance. We conclude by showing how some non-independence Metropolis–Hastings algorithms can perform better than the IMS and deriving general lower and upper bounds for the mean first hitting times of a Metropolis–Hastings algorithm.     

Convergence analysis of the projective scaling algorithm based on a long-step homogeneous affine scaling algorithm

In this paper we give a new convergence analysis of a projective scaling algorithm. We consider a long-step affine scaling algorithm applied to a homogeneous linear programming problem obtained from the original linear programming problem. This algorithm takes a fixed fraction λ≤2/3 of the way towards the boundary of the nonnegative orthant at each iteration. The iteration sequence for the original problem is obtained by pulling back the homogeneous iterates onto the original feasible region with a conical projection, which generates the same search direction as the original projective scaling algorithm at each iterate. The recent convergence results for the long-step affine scaling algorithm by the authors are applied to this algorithm to obtain some convergence results on the projective scaling algorithm. Specifically, we will show (i) polynomiality of the algorithm with complexities of O(nL) and O(n ² L) iterations for λ<2/3 and λ=2/3, respectively; (ii) global covnergence of the algorithm when the optimal face is unbounded; (iii) convergence of the primal iterates to a relative interior point of the optimal face; (iv) convergence of the dual estimates to the analytic center of the dual optimal face; and (v) convergence of the reduction rate of the objective function value to 1−λ.     

A semiparametric density estimator based on elliptical distributions

In the paper we study a semiparametric density estimation method based on the model of an elliptical distribution. The method considered here shows a way to overcome problems arising from the curse of dimensionality. The optimal rate of the uniform strong convergence of the estimator under consideration coincides with the optimal rate for the usual one-dimensional kernel density estimator except in a neighbourhood of the mean. Therefore the optimal rate does not depend on the dimension. Moreover, asymptotic normality of the estimator is proved.