Markov Chain Sampling Methods for Dirichlet Process Mixture Models

Abstract

This article reviews Markov chain methods for sampling from the posterior distribution of a Dirichlet process mixture model and presents two new classes of methods. One new approach is to make Metropolis—Hastings updates of the indicators specifying which mixture component is associated with each observation, perhaps supplemented with a partial form of Gibbs sampling. The other new approach extends Gibbs sampling for these indicators by using a set of auxiliary parameters. These methods are simple to implement and are more efficient than previous ways of handling general Dirichlet process mixture models with non-conjugate priors.     

Statistical analysis of hierarchical stochastic models: Examples and approaches

This paper introduces and illustrates the concept of hierarchical or random parameter stochastic process models. These models arise when members of a population each generate a stochastic process governed by certain parameters and the values of the parameters may be viewed as single realizations of random variables. The paper treats the estimation of the individual parameter values and the parameters of the superpopulation distribution. Examples from system reliability, pharmacokinetic compartment models, and criminal careers are introduced; a reliability (Poisson process-exponential interval) process is examined in greater detail. An explicit, approximate, robust estimator of individual (log) failure rates is presented for the case of a long-tailed (Studentt) superpopulation. This estimator exhibits desirable limited shrinkage properties, refusing to borrow unjustified strength. Numerical properties of such estimators are described more fully elsewhere.     

Laplace Expansions in Markov chain Monte Carlo Algorithms

Complex hierarchical models lead to a complicated likelihood and then, in a Bayesian analysis, to complicated posterior distributions. To obtain Bayes estimates such as the posterior mean or Bayesian confidence regions, it is therefore necessary to simulate the posterior distribution using a method such as an MCMC algorithm. These algorithms often get slower as the number of observations increases, especially when the latent variables are considered. To improve the convergence of the algorithm, we propose to decrease the number of parameters to simulate at each iteration by using a Laplace approximation on the nuisance parameters. We provide a theoretical study of the impact that such an approximation has on the target posterior distribution. We prove that the distance between the true target distribution and the approximation becomes of order O(N^?a) with a ∈ (0, 1), a close to 1, as the number of observations N increases. A simulation study illustrates the theoretical results. The approximated MCMC algorithm behaves extremely well on an example which is driven by a study on HIV patients.     

Functional Compatibility,Markov Chains,and Gibbs Sampling with Improper Posteriors

Abstract

The members of a set of conditional probability density functions are called compatible if there exists a joint probability density function that generates them. We generalize this concept by calling the conditionals functionally compatible if there exists a non-negative function that behaves like a joint density as far as generating the conditionals according to the probability calculus, but whose integral over the whole space is not necessarily finite. A necessary and sufficient condition for functional compatibility is given that provides a method of calculating this function, if it exists. A Markov transition function is then constructed using a set of functionally compatible conditional densities and it is shown, using the compatibility results, that the associated Markov chain is positive recurrent if and only if the conditionals are compatible. A Gibbs Markov chain, constructed via “Gibbs conditionals” from a hierarchical model with an improper posterior, is a special case. Therefore, the results of this article can be used to evaluate the consequences of applying the Gibbs sampler when the posterior's impropriety is unknown to the user. Our results cannot, however, be used to detect improper posteriors. Monte Carlo approximations based on Gibbs chains are shown to have undesirable limiting behavior when the posterior is improper. The results are applied to a Bayesian hierarchical one-way random effects model with an improper posterior distribution. The model is simple, but also quite similar to some models with improper posteriors that have been used in conjunction with the Gibbs sampler in the literature.     

Approximated Bayes and empirical Bayes confidence intervals—The known variance case

In this paper hierarchical Bayes and empirical Bayes results are used to obtain confidence intervals of the population means in the case of real problems. This is achieved by approximating the posterior distribution with a Pearson distribution. In the first example hierarchical Bayes confidence intervals for the Efron and Morris (1975, J. Amer. Statist. Assoc., 70, 311–319) baseball data are obtained. The same methods are used in the second example to obtain confidence intervals of treatment effects as well as the difference between treatment effects in an analysis of variance experiment. In the third example hierarchical Bayes intervals of treatment effects are obtained and compared with normal approximations in the unequal variance case.Financially supported by the CSIR and the University of the Orange Free State, Central Research Fund.     

Maximum Likelihood for Generalized Linear Models with Nested Random Effects via High-Order,Multivariate Laplace Approximation

Abstract

Nested random effects models are often used to represent similar processes occurring in each of many clusters. Suppose that, given cluster-specific random effects b, the data y are distributed according to f(y|b, Θ), while b follows a density p(b|Θ). Likelihood inference requires maximization of ∫ f(y|b, Θ)p(b|Θdb with respect to Θ. Evaluation of this integral often proves difficult, making likelihood inference difficult to obtain. We propose a multivariate Taylor series approximation of the log of the integrand that can be made as accurate as desired if the integrand and all its partial derivatives with respect to b are continuous in the neighborhood of the posterior mode of b|Θ,y. We then apply a Laplace approximation to the integral and maximize the approximate integrated likelihood via Fisher scoring. We develop computational formulas that implement this approach for two-level generalized linear models with canonical link and multivariate normal random effects. A comparison with approximations based on penalized quasi-likelihood, Gauss—Hermite quadrature, and adaptive Gauss-Hermite quadrature reveals that, for the hierarchical logistic regression model under the simulated conditions, the sixth-order Laplace approach is remarkably accurate and computationally fast.     

Reparameterized and Marginalized Posterior and Predictive Sampling for Complex Bayesian Geostatistical Models

This article proposes a four-pronged approach to efficient Bayesian estimation and prediction for complex Bayesian hierarchical Gaussian models for spatial and spatiotemporal data. The method involves reparameterizing the covariance structure of the model, reformulating the means structure, marginalizing the joint posterior distribution, and applying a simplex-based slice sampling algorithm. The approach permits fusion of point-source data and areal data measured at different resolutions and accommodates nonspatial correlation and variance heterogeneity as well as spatial and/or temporal correlation. The method produces Markov chain Monte Carlo samplers with low autocorrelation in the output, so that fewer iterations are needed for Bayesian inference than would be the case with other sampling algorithms. Supplemental materials are available online.     

Adaptive Markov Chain Monte Carlo for Bayesian Variable Selection

We describe adaptive Markov chain Monte Carlo (MCMC) methods for sampling posterior distributions arising from Bayesian variable selection problems. Point-mass mixture priors are commonly used in Bayesian variable selection problems in regression. However, for generalized linear and nonlinear models where the conditional densities cannot be obtained directly, the resulting mixture posterior may be difficult to sample using standard MCMC methods due to multimodality. We introduce an adaptive MCMC scheme that automatically tunes the parameters of a family of mixture proposal distributions during simulation. The resulting chain adapts to sample efficiently from multimodal target distributions. For variable selection problems point-mass components are included in the mixture, and the associated weights adapt to approximate marginal posterior variable inclusion probabilities, while the remaining components approximate the posterior over nonzero values. The resulting sampler transitions efficiently between models, performing parameter estimation and variable selection simultaneously. Ergodicity and convergence are guaranteed by limiting the adaptation based on recent theoretical results. The algorithm is demonstrated on a logistic regression model, a sparse kernel regression, and a random field model from statistical biophysics; in each case the adaptive algorithm dramatically outperforms traditional MH algorithms. Supplementary materials for this article are available online.     

Method of Moments Using Monte Carlo Simulation

Abstract

We present a computational approach to the method of moments using Monte Carlo simulation. Simple algebraic identities are used so that all computations can be performed directly using simulation draws and computation of the derivative of the log-likelihood. We present a simple implementation using the Newton-Raphson algorithm with the understanding that other optimization methods may be used in more complicated problems. The method can be applied to families of distributions with unknown normalizing constants and can be extended to least squares fitting in the case that the number of moments observed exceeds the number of parameters in the model. The method can be further generalized to allow “moments” that are any function of data and parameters, including as a special case maximum likelihood for models with unknown normalizing constants or missing data. In addition to being used for estimation, our method may be useful for setting the parameters of a Bayes prior distribution by specifying moments of a distribution using prior information. We present two examples—specification of a multivariate prior distribution in a constrained-parameter family and estimation of parameters in an image model. The former example, used for an application in pharmacokinetics, motivated this work. This work is similar to Ruppert's method in stochastic approximation, combines Monte Carlo simulation and the Newton-Raphson algorithm as in Penttinen, uses computational ideas and importance sampling identities of Gelfand and Carlin, Geyer, and Geyer and Thompson developed for Monte Carlo maximum likelihood, and has some similarities to the maximum likelihood methods of Wei and Tanner.     

Bayesian joint quantile regression for mixed effects models with censoring and errors in covariates

In this paper, we discuss Bayesian joint quantile regression of mixed effects models with censored responses and errors in covariates simultaneously using Markov Chain Monte Carlo method. Under the assumption of asymmetric Laplace error distribution, we establish a Bayesian hierarchical model and derive the posterior distributions of all unknown parameters based on Gibbs sampling algorithm. Three cases including multivariate normal distribution and other two heavy-tailed distributions are considered for fitting random effects of the mixed effects models. Finally, some Monte Carlo simulations are performed and the proposed procedure is illustrated by analyzing a group of AIDS clinical data set.     

Asymptotic Corrections for Multivariate Posterior Moments with Factored Likelihood Functions

Abstract

Asymptotic corrections are used to compute the means and the variance-covariance matrix of multivariate posterior distributions that are formed from a normal prior distribution and a likelihood function that factors into separate functions for each variable in the posterior distribution. The approximations are illustrated using data from the National Assessment of Educational Progress (NAEP). These corrections produce much more accurate approximations than those produced by two different normal approximations. In a second potential application, the computational methods are applied to logistic regression models for severity adjustment of hospital-specific mortality rates.     

A Comparison of Two Sequential Metropolis-Hastings Algorithms with Standard Simulation Techniques in Bayesian Inference in Reliability Models Involving a Generalized Gamma Distribution

In this paper we consider the generalized gamma distribution as introduced in Gåsemyr and Natvig (1998). This distribution enters naturally in Bayesian inference in exponential survival models with left censoring. In the paper mentioned above it is shown that the weighted sum of products of generalized gamma distributions is a conjugate prior for the parameters of component lifetimes, having autopsy data in a Marshall-Olkin shock model. A corresponding result is shown in Gåsemyr and Natvig (1999) for independent, exponentially distributed component lifetimes in a model with partial monitoring of components with applications to preventive system maintenance. A discussion in the present paper strongly indicates that expressing the posterior distribution in terms of the generalized gamma distribution is computationally efficient compared to using the ordinary gamma distribution in such models. Furthermore, we present two types of sequential Metropolis-Hastings algorithms that may be used in Bayesian inference in situations where exact methods are intractable. Finally these types of algorithms are compared with standard simulation techniques and analytical results in arriving at the posterior distribution of the parameters of component lifetimes in special cases of the mentioned models. It seems that one of these types of algorithms may be very favorable when prior assessments are updated by several data sets and when there are significant discrepancies between the prior assessments and the data.     

A review of multi-product pricing models

The importance of good pricing strategies in business theory is clearly recognized, as can be seen from the huge volume of pricing research done over the years. What we attempt to do is to provide a general review of multi-product pricing models, focusing primarily on those where demands are explicitly dependent on prices. As the pricing decision may be made jointly with other economic parameters, we will not only review models that focus solely on pricing; we will also discuss models where pricing choices are made jointly with other decisions like production or distribution of resources.     

Bayesian inference in spherical linear models: robustness and conjugate analysis

The early work of Zellner on the multivariate Student-t linear model has been extended to Bayesian inference for linear models with dependent non-normal error terms, particularly through various papers by Osiewalski, Steel and coworkers. This article provides a full Bayesian analysis for a spherical linear model. The density generator of the spherical distribution is here allowed to depend both on the precision parameter φ and on the regression coefficients β. Another distinctive aspect of this paper is that proper priors for the precision parameter are discussed.The normal-chi-squared family of prior distributions is extended to a new class, which allows the posterior analysis to be carried out analytically. On the other hand, a direct joint modelling of the data vector and of the parameters leads to conjugate distributions for the regression and the precision parameters, both individually and jointly. It is shown that some model specifications lead to Bayes estimators that do not depend on the choice of the density generator, in agreement with previous results obtained in the literature under different assumptions. Finally, the distribution theory developed to tackle the main problem is useful on its own right.     

Estimating Functions for Nonlinear Time Series Models

This paper discusses the problem of estimation for two classes of nonlinear models, namely random coefficient autoregressive (RCA) and autoregressive conditional heteroskedasticity (ARCH) models. For the RCA model, first assuming that the nuisance parameters are known we construct an estimator for parameters of interest based on Godambe's asymptotically optimal estimating function. Then, using the conditional least squares (CLS) estimator given by Tjøstheim (1986, Stochastic Process. Appl., 21, 251–273) and classical moment estimators for the nuisance parameters, we propose an estimated version of this estimator. These results are extended to the case of vector parameter. Next, we turn to discuss the problem of estimating the ARCH model with unknown parameter vector. We construct an estimator for parameters of interest based on Godambe's optimal estimator allowing that a part of the estimator depends on unknown parameters. Then, substituting the CLS estimators for the unknown parameters, the estimated version is proposed. Comparisons between the CLS and estimated optimal estimator of the RCA model and between the CLS and estimated version of the ARCH model are given via simulation studies.     

On the Convergence of Successive Substitution Sampling

Abstract

The problem of finding marginal distributions of multidimensional random quantities has many applications in probability and statistics. Many of the solutions currently in use are very computationally intensive. For example, in a Bayesian inference problem with a hierarchical prior distribution, one is often driven to multidimensional numerical integration to obtain marginal posterior distributions of the model parameters of interest. Recently, however, a group of Monte Carlo integration techniques that fall under the general banner of successive substitution sampling (SSS) have proven to be powerful tools for obtaining approximate answers in a very wide variety of Bayesian modeling situations. Answers may also be obtained at low cost, both in terms of computer power and user sophistication. Important special cases of SSS include the “Gibbs sampler” described by Gelfand and Smith and the “IP algorithm” described by Tanner and Wong. The major problem plaguing users of SSS is the difficulty in ascertaining when “convergence” of the algorithm has been obtained. This problem is compounded by the fact that what is produced by the sampler is not the functional form of the desired marginal posterior distribution, but a random sample from this distribution. This article gives a general proof of the convergence of SSS and the sufficient conditions for both strong and weak convergence, as well as a convergence rate. We explore the connection between higher-order eigenfunctions of the transition operator and accelerated convergence via good initial distributions. We also provide asymptotic results for the sampling component of the error in estimating the distributions of interest. Finally, we give two detailed examples from familiar exponential family settings to illustrate the theory.     

Bayesian‐type count data models with varying coefficients: estimation and testing in the presence of overdispersion

In this paper we study varying‐coefficient models for count data. A Bayesian approach is taken to model the variability of the regression parameters. Based on a Kalman filter procedure the varying coefficients are estimated as the mode of the posterior distribution. All hyperparameters, including an overdispersion parameter in the negative binomial varying‐coefficient model (NBVC), are estimated as ML‐estimators using an EM‐type algorithm. A bootstrapping test of the fixed‐coefficient hypothesis against a varying‐coefficient alternative is proposed, which is evaluated running a simulation study. The study shows that the choice of a suitable count data model is of special importance in the framework of varying‐coefficient models. The methodology is illustrated analysing the determinants of the number of individual doctor visits. Copyright © 2001 John Wiley & Sons, Ltd.     

Multivariate spatial regression models

This paper describes the inference procedures required to perform Bayesian inference to some multivariate econometric models. These models have a spatial component built into commonly used multivariate models. In particular, the common component models are addressed and extended to accommodate for spatial dependence. Inference procedures are based on a variety of simulation-based schemes designed to obtain samples from the posterior distribution of model parameters. They are also used to provide a basis to forecast new observations.     

Hierarchical linear regression models for conditional quantiles

The quantile regression has several useful features and therefore is gradually developing into a comprehensive approach to the statistical analysis of linear and nonlinear response models, but it cannot deal effectively with the data with a hierarchical structure. In practice, the existence of such data hierarchies is neither accidental nor ignorable, it is a common phenomenon. To ignore this hierarchical data structure risks overlooking the importance of group effects, and may also render many of the traditional statistical analysis techniques used for studying data relationships invalid. On the other hand, the hierarchical models take a hierarchical data structure into account and have also many applications in statistics, ranging from overdispersion to constructing min-max estimators. However, the hierarchical models are virtually the mean regression, therefore, they cannot be used to characterize the entire conditional distribution of a dependent variable given high-dimensional covariates. Furthermore, the estimated coefficient vector (marginal effects) is sensitive to an outlier observation on the dependent variable. In this article, a new approach, which is based on the Gauss-Seidel iteration and taking a full advantage of the quantile regression and hierarchical models, is developed. On the theoretical front, we also consider the asymptotic properties of the new method, obtaining the simple conditions for an n1/2-convergence and an asymptotic normality. We also illustrate the use of the technique with the real educational data which is hierarchical and how the results can be explained.