Sampling Correlation Matrices in Bayesian Models With Correlated Latent Variables

Hierarchical model specifications using latent variables are frequently used to reflect correlation structure in data. Motivated by the structure of a Bayesian multivariate probit model, we demonstrate a parameter-extended Metropolis-Hastings algorithm for sampling from the posterior distribution of a correlation matrix. Our sampling algorithms lead directly to two readily interpretable families of prior distributions for a correlation matrix. The methodology is illustrated through a simulation study and through an application with repeated binary outcomes on individuals from a study of a suicide prevention intervention.     

High Breakdown Estimation for Multiple Populations with Applications to Discriminant Analysis

We consider S-estimators of multivariate location and common dispersion matrix in multiple populations. Instead of averaging the robust estimates of the individual covariance matrices, as used by Todorov, Neykov and Neytchev (1990), the observations are pooled for estimating the common covariance more efficiently. Two such proposals are evaluated by a breakdown point analysis and Monte Carlo simulations. Their applications to the discriminant analysis are also considered.     

Efficient algorithms for generating truncated multivariate normal distributions

Sampling from a truncated multivariate normal distribution (TMVND) constitutes the core computational module in fitting many statistical and econometric models. We propose two efficient methods, an iterative data augmentation (DA) algorithm and a non-iterative inverse Bayes formulae (IBF) sampler, to simulate TMVND and generalize them to multivariate normal distributions with linear inequality constraints. By creating a Bayesian incomplete-data structure, the posterior step of the DA algorithm directly generates random vector draws as opposed to single element draws, resulting obvious computational advantage and easy coding with common statistical software packages such as S-PLUS, MATLAB and GAUSS. Furthermore, the DA provides a ready structure for implementing a fast EM algorithm to identify the mode of TMVND, which has many potential applications in statistical inference of constrained parameter problems. In addition, utilizing this mode as an intermediate result, the IBF sampling provides a novel alternative to Gibbs sampling and eliminates problems with convergence and possible slow convergence due to the high correlation between components of a TMVND. The DA algorithm is applied to a linear regression model with constrained parameters and is illustrated with a published data set. Numerical comparisons show that the proposed DA algorithm and IBF sampler are more efficient than the Gibbs sampler and the accept-reject algorithm.     

The affine equivariant sign covariance matrix: asymptotic behavior and efficiencies

We consider the affine equivariant sign covariance matrix (SCM) introduced by Visuri et al. (J. Statist. Plann. Inference 91 (2000) 557). The population SCM is shown to be proportional to the inverse of the regular covariance matrix. The eigenvectors and standardized eigenvalues of the covariance matrix can thus be derived from the SCM. We also construct an estimate of the covariance and correlation matrix based on the SCM. The influence functions and limiting distributions of the SCM and its eigenvectors and eigenvalues are found. Limiting efficiencies are given in multivariate normal and t-distribution cases. The estimates are highly efficient in the multivariate normal case and perform better than estimates based on the sample covariance matrix for heavy-tailed distributions. Simulations confirmed these findings for finite-sample efficiencies.     

Efficient generation of random nonsingular matrices

We present an efficient algorithm for generating an n × n nonsingular matrix uniformly over a finite field. This algorithm is useful for several cryptographic and checking applications. Over GF[2] our algorithm runs in expected time M(n) + O(n²), where M(n) is the time needed to multiply two n × n matrices, and the expected number of random bits it uses is n² + 3. (Over other finite fields we use n² + O(1) random field elements on average.) This is more efficient than the standard method for solving this problem, both in terms of expected running time and the expected number of random bits used. The standard method is to generate random n × n matrices until we produce one with nonzero determinant. In contrast, our technique directly produces a random matrix guaranteed to have nonzero determinant. We also introduce efficient algorithms for related problems such as uniformly generating singular matrices or matrices with fixed determinant. © 1993 John Wiley & Sons, Inc.     

MCMC Estimation of Restricted Covariance Matrices

This article is motivated by the difficulty of applying standard simulation techniques when identification constraints or theoretical considerations induce covariance restrictions in multivariate models. To deal with this difficulty, we build upon a decomposition of positive definite matrices and show that it leads to straightforward Markov chain Monte Carlo samplers for restricted covariance matrices. We introduce the approach by reviewing results for multivariate Gaussian models without restrictions, where standard conjugate priors on the elements of the decomposition induce the usual Wishart distribution on the precision matrix and vice versa. The unrestricted case provides guidance for constructing efficient Metropolis–Hastings and accept-reject Metropolis–Hastings samplers in more complex settings, and we describe in detail how simulation can be performed under several important constraints. The proposed approach is illustrated in a simulation study and two applications in economics. Supplemental materials for this article (appendixes, data, and computer code) are available online.     

Efficient Covariance Approximations for Large Sparse Precision Matrices

The use of sparse precision (inverse covariance) matrices has become popular because they allow for efficient algorithms for joint inference in high-dimensional models. Many applications require the computation of certain elements of the covariance matrix, such as the marginal variances, which may be nontrivial to obtain when the dimension is large. This article introduces a fast Rao–Blackwellized Monte Carlo sampling-based method for efficiently approximating selected elements of the covariance matrix. The variance and confidence bounds of the approximations can be precisely estimated without additional computational costs. Furthermore, a method that iterates over subdomains is introduced, and is shown to additionally reduce the approximation errors to practically negligible levels in an application on functional magnetic resonance imaging data. Both methods have low memory requirements, which is typically the bottleneck for competing direct methods.     

Efficient Bayesian Inference for Multivariate Probit Models With Sparse Inverse Correlation Matrices

We propose a Bayesian approach for inference in the multivariate probit model, taking into account the association structure between binary observations. We model the association through the correlation matrix of the latent Gaussian variables. Conditional independence is imposed by setting some off-diagonal elements of the inverse correlation matrix to zero and this sparsity structure is modeled using a decomposable graphical model. We propose an efficient Markov chain Monte Carlo algorithm relying on a parameter expansion scheme to sample from the resulting posterior distribution. This algorithm updates the correlation matrix within a simple Gibbs sampling framework and allows us to infer the correlation structure from the data, generalizing methods used for inference in decomposable Gaussian graphical models to multivariate binary observations. We demonstrate the performance of this model and of the Markov chain Monte Carlo algorithm on simulated and real datasets. This article has online supplementary materials.     

Efficient strategies for deriving the subset VAR models

Algorithms for computing the subset Vector Autoregressive (VAR) models are proposed. These algorithms can be used to choose a subset of the most statistically-significant variables of a VAR model. In such cases, the selection criteria are based on the residual sum of squares or the estimated residual covariance matrix. The VAR model with zero coefficient restrictions is formulated as a Seemingly Unrelated Regressions (SUR) model. Furthermore, the SUR model is transformed into one of smaller size, where the exogenous matrices comprise columns of a triangular matrix. Efficient algorithms which exploit the common columns of the exogenous matrices, sparse structure of the variance-covariance of the disturbances and special properties of the SUR models are investigated. The main computational tool of the selection strategies is the generalized QR decomposition and its modification. This work is in part supported by the Swiss National Foundation for Research Grants 101312-100757/1, 200020-10016/1 and 101412-105978. Correspondence to: Cristian Gatu, Institut d'informatique, Université de Neuchatel, Emile-Argand 11, Case Postale 2, CH-2007 Neuchatel, Switzerland. e-mail: erricos@ucy.ac.cy     

Computing powers of arbitrary Hessenberg matrices

An efficient algorithm for the computation of powers of an n × n arbitrary lower Hessenberg matrix is presented. Numerical examples are used to show the computational details. A comparison of the algorithm with two other methods of matrix multiplication proposed by Brent and by Winograd is included. Related algorithms were proposed earlier by Datta and Datta for lower Hessenberg matrices with unit super-diagonal elements, and by Barnett for the companion matrix.     

Testing some covariance structures under a growth curve model

This paper considers three types of problems: (i) the problem of independence of two sets, (ii) the problem of sphericity of the covariance matrix Σ, and (iii) the problem of intraclass model for the covariance matrix Σ, when the column vectors of X are independently distributed as multivariate normal with covariance matrix Σ and E(X) = BξA,A and B being given matrices and ξ and Σ being unknown. These problems are solved by the likelihood ratio test procedures under some restrictions on the models, and the null distributions of the test statistics are established.     

Covariance Decompositions for Accurate Computation in Bayesian Scale-Usage Models

Analyses of multivariate ordinal probit models typically use data augmentation to link the observed (discrete) data to latent (continuous) data via a censoring mechanism defined by a collection of “cutpoints.” Most standard models, for which effective Markov chain Monte Carlo (MCMC) sampling algorithms have been developed, use a separate (and independent) set of cutpoints for each element of the multivariate response. Motivated by the analysis of ratings data, we describe a particular class of multivariate ordinal probit models where it is desirable to use a common set of cutpoints. While this approach is attractive from a data-analytic perspective, we show that the existing efficient MCMC algorithms can no longer be accurately applied. Moreover, we show that attempts to implement these algorithms by numerically approximating required multivariate normal integrals over high-dimensional rectangular regions can result in severely degraded estimates of the posterior distribution. We propose a new data augmentation that is based on a covariance decomposition and that admits a simple and accurate MCMC algorithm. Our data augmentation requires only that univariate normal integrals be evaluated, which can be done quickly and with high accuracy. We provide theoretical results that suggest optimal decompositions within this class of data augmentations, and, based on the theory, recommend default decompositions that we demonstrate work well in practice. This article has supplementary material online.     

A note on the <Emphasis Type="BoldItalic">O</Emphasis>(<Emphasis Type="BoldItalic">n</Emphasis>)-storage implementation of the GKO algorithm and its adaptation to Trummer-like matrices

We propose a new O(n)-space implementation of the GKO-Cauchy algorithm for the solution of linear systems where the coefficient matrix is Cauchy-like. Moreover, this new algorithm makes a more efficient use of the processor cache memory; for matrices of size larger than n ≈ 500–1,000, it outperforms the customary GKO algorithm. We present an applicative case of Cauchy-like matrices with non-reconstructible main diagonal. In this special instance, the O(n) space algorithms can be adapted nicely to provide an efficient implementation of basic linear algebra operations in terms of the low displacement-rank generators.     

Approximations to the distribution of the sample correlation matrix

In this article, multivariate density expansions for the sample correlation matrix R are derived. The density of R is expressed through multivariate normal and through Wishart distributions. Also, an asymptotic expansion of the characteristic function of R is derived and the main terms of the first three cumulants of R are obtained in matrix form. These results make it possible to obtain asymptotic density expansions of multivariate functions of R in a direct way.     

An Interior-Point Method for Approximate Positive Semidefinite Completions

Given a nonnegative, symmetric matrix of weights, H, we study the problem of finding an Hermitian, positive semidefinite matrix which is closest to a given Hermitian matrix, A, with respect to the weighting H. This extends the notion of exact matrix completion problems in that, H _ij =0 corresponds to the element A _ij being unspecified (free), while H _ij large in absolute value corresponds to the element A _ij being approximately specified (fixed).We present optimality conditions, duality theory, and two primal-dual interior-point algorithms. Because of sparsity considerations, the dual-step-first algorithm is more efficient for a large number of free elements, while the primal-step-first algorithm is more efficient for a large number of fixed elements.Included are numerical tests that illustrate the efficiency and robustness of the algorithms     

Generating multivariate correlated samples

Summary  Generating random samples from multivariate distributions is a common, requirement in many fields of study. Often the complete joint distribution is not specified to the scientist. This paper addresses the situation where only the marginals and the correlation matrix are specified. We suggest a deterministic algorithm, PERMCORR, to approximately achieve the required correlation structure that can be used to get good initial values to standard stochastic algorithms. In many situations the output of PERMCORR is already accurate enough to preempt any need for running an expensive stochastic algorithm. We provide some theoretical justification for our method as well as simulation studies. We also provide a bootstrap technique based on PERMCORR.     

Bayesian Prediction Analysis for Growth Curve Model Using Noninformative Priors

We apply a Bayesian approach to the problem of prediction in an unbalanced growth curve model using noninformative priors. Due to the complexity of the model, no analytic forms of the predictive densities are available. We propose both approximations and a prediction-oriented Metropolis-Hastings sampling algorithm for two types of prediction, namely the prediction of future observations for a new subject and the prediction of future values for a partially observed subject. They are illustrated and compared through real data and simulation studies. Two of the approximations compare favorably with the approximation in Fearn (1975, Biometrika, 62, 89–100) and are very comparable to the more accurate Rao-Blackwellization from Metropolis-Hastings sampling algorithm.     

Fast block Toeplitz orthogonalization

Summary Algorithms are presented which compute theQR factorization of a block-Toeplitz matrix inO(n) ² block-operations, wheren is the block-order of the matrix and a block-operation is a multiplication, inversion or a set of Householder operations involving one or two blocks. The algorithms are in general analogous to those presented in the scalar Toeplitz case in a previous paper, but the basic operation is the Householder transform rather than the Givens transform, and the computation of the Householder coefficients and other working matrices requires special treatment. Two algorithms are presented-the first computes onlyR explicitly and the second computes bothQ andR.     

An alternating direction algorithm for matrix completion with nonnegative factors

This paper introduces an algorithm for the nonnegative matrix factorization-and-completion problem, which aims to find nonnegative low-rank matrices X and Y so that the product XY approximates a nonnegative data matrix M whose elements are partially known (to a certain accuracy). This problem aggregates two existing problems: (i) nonnegative matrix factorization where all entries of M are given, and (ii) low-rank matrix completion where nonnegativity is not required. By taking the advantages of both nonnegativity and low-rankness, one can generally obtain superior results than those of just using one of the two properties. We propose to solve the non-convex constrained least-squares problem using an algorithm based on the classical alternating direction augmented Lagrangian method. Preliminary convergence properties of the algorithm and numerical simulation results are presented. Compared to a recent algorithm for nonnegative matrix factorization, the proposed algorithm produces factorizations of similar quality using only about half of the matrix entries. On tasks of recovering incomplete grayscale and hyperspectral images, the proposed algorithm yields overall better qualities than those produced by two recent matrix-completion algorithms that do not exploit nonnegativity.