We present methods for the computation of multivariate normal probabilities on parallel/ distributed systems. After a transformation of the initial integral, an approximation can be obtained using Monte-Carlo or quasirandom methods. We propose a meta-algorithm for asynchronous sampling methods and derive efficient parallel algorithms for the computation of MVN distribution functions, including a method based on randomized Korobov and Richtmyer sequences. Timing results of the implementations using the MPI parallel environment are given. 1 Introduction The computation of the multivariate normal distribution function F (a; b) = j\Sigmaj \Gamma 1 2 (2) \Gamma n 2 Z b a e \Gamma 1 2 x \Sigma \Gamma1 x dx: (1) often leads to computational-intensive integration problems. Here \Sigma is an n \Theta n symmetric positive definite covariance matrix; furthermore one of the limits in each integration variable may be infinite. Genz [5] performs a sequence of transformations resu...

This paper provides a practical simulation-based Bayesian and non-Bayesian analysis of correlated binary data using the multivariate probit model. The posterior distribution is simulated by Markov chain Monte Carlo methods and maximum likelihood estimates are obtained by a Monte Carlo version of the EM algorithm. A practical approach for the computation of Bayes factors from the simulation output is also developed. The methods are applied to a dataset with a bivariate binary response, to a four-year longitudinal dataset from the Six Cities study of the health effects of air pollution and to a sevenvariate binary response dataset on the labour supply of married women from the Panel Survey of Income Dynamics.

This paper compares methods for the numerical computation of multivariate t-probabilities for hyperrectangular integration regions. Methods based on acceptance-rejection, spherical-radial transformations and separation-of-variables transformations are considered. Tests using randomly chosen problems show that the most efficient numerical methods use a transformation developed by Genz (1992) for multivariate normal probabilities. These methods allow moderately accurate multivariate t-probabilities to be quickly computed for problems with as many as twenty variables. Methods for the non-central multivariate t-distribution are also described. Key Words: multivariate t-distribution, non-central distribution, numerical integration, statistical computation. 1 Introduction A common problem in many statistics applications is the numerical computation of the multivariate t (MVT) distribution function (see Tong, 1990) defined by T(a; b; \Sigma; ) = \Gamma( +m 2 ) \Gamma( 2 ) p j\Sigma...

Integration lattices are one of the main types of low discrepancy sets used in quasi-Monte Carlo methods. However, they have the disadvantage of being of fixed size. This article describes the construction of an infinite sequence of points, the first b m of which form a lattice for any non-negative integer m. Thus, if the quadrature error using an initial lattice is too large, the lattice can be extended without discarding the original points. Generating vectors for extensible lattices are found by minimizing a loss function based on some measure of discrepancy or nonuniformity of the lattice. The spectral test used for finding pseudo-random number generators is one important example of such a discrepancy. The performance of the extensible lattices proposed here is compared to that of other methods for some practical quadrature problems.

. Correlation matrices---symmetric positive semidefinite matrices with unit diagonal--- are important in statistics and in numerical linear algebra. For simulation and testing it is desirable to be able to generate random correlation matrices with specified eigenvalues (which must be nonnegative and sum to the dimension of the matrix). A popular algorithm of Bendel and Mickey takes a matrix having the specified eigenvalues and uses a finite sequence of Given rotations to introduce 1s on the diagonal. We give improved formulae for computing the rotations and prove that the resulting algorithm is numerically stable. We show by example that the formulae originally proposed, which are used in certain existing Fortran implementations, can lead to serious instability. We also show how to modify the algorithm to generate a rectangular matrix with columns of unit 2-norm. Such a matrix represents a correlation matrix in factored form, which can be preferable to representing the matrix itself, ...

Models of dynamical systems based on predictive state representations (PSRs) are defined strictly in terms of observable quantities, in contrast with traditional models (such as Hidden Markov Models) that use latent variables or statespace representations. In addition, PSRs have an effectively infinite memory, allowing them to model some systems that finite memory-based models cannot. Thus far, PSR models have primarily been developed for domains with discrete observations. Here, we develop the Predictive Linear-Gaussian (PLG) model, a class of PSR models for domains with continuous observations. We show that PLG models subsume Linear Dynamical System models (also called Kalman filter models or state-space models) while using fewer parameters. We also introduce an algorithm to estimate PLG parameters from data, and contrast it with standard Expectation Maximization (EM) algorithms used to estimate Kalman filter parameters. We show that our algorithm is a consistent estimation procedure and present preliminary empirical results suggesting that our algorithm outperforms EM, particularly as the model dimension increases. 1

We consider the problem of comparing a finite number of stochastic systems with respect to a single system (designated as the "standard") via simulation experiments. The comparison is based on expected performance, and the goal is to determine if any system has larger expected performance than the standard, and if so to identify the best of the alternatives. In this paper we provide two-stage experiment design and analysis procedures to solve the problem for a variety of scenarios, including when we encounter unequal variances across systems, and when we use the variance reduction technique of common random numbers and it is appropriate to do so. The emphasis is added because in some cases common random numbers can be counterproductive when performing comparisons with a standard. We also provide methods for estimating the critical constants required by our procedures, present a portion of an extensive empirical study and demonstrate one of the procedures via a numerical example. 1 Intr...

. Multivariate normal probabilities, which are used for statistical inference, must be computed numerically. This article describes a new rank-1 lattice quadrature rule and its application to computing multivariate normal probabilities. In contrast to existing lattice rules the number of integrand evaluations need not be specified in advance. When compared to existing algorithms for computing multivariate normal probabilities the new algorithm is more efficient when high accuracy is required and/or the number of variables is large. 1 Introduction The most important probability distribution is the Gaussian or normal probability distribution. Normal probabilities are used to perform statistical inference and construct confidence intervals. The definition of the normal probability distribution involves an integral which cannot be evaluated in terms of elementary functions. Therefore, numerical methods are needed. Many software packages contain routines for evaluating univariate normal pr...

In this paper, we present new algorithms that can replace the diagonal entries of a Hermitian matrix by any set of diagonal entries that majorize the original set without altering the eigenvalues of the matrix. They perform this feat by applying a sequence of (N − 1) or fewer plane rotations, where N is the dimension of the matrix. Both the Bendel–Mickey and the Chan–Li algorithms are special cases of the proposed procedures. Using the fact that a positive semidefinite matrix can always be factored as X ∗ X, we also provide more efficient versions of the algorithms that can directly construct factors with specified singular values and column norms. We conclude with some open problems related to the construction of Hermitian matrices with joint diagonal and spectral properties.

In this paper, we use multivariate logistic regression models to incorporate correlation among binary response data. Our objective is to develop a variable subset selection procedure to identify important covariates in predicting correlated binary responses using a Bayesian approach. In order to incorporate available prior information, we propose a class of informative prior distributions on the model parameters and on the model space. The propriety of the proposed informative prior is investigated in detail. Novel computational algorithms are also developed for sampling from the posterior distribution as well as for computing posterior model probabilities. Finally, a simulated data example and a real data example from a prostate cancer study are used to illustrate the proposed methodology.