# [PERSON] and Mendell-Elston Estimation of the High-Dimensional Multivariate Normal Distribution
## 1 Introduction  
In applied multivariate statistical analysis one is frequently faced with the problem of estimating the multivariate normal (MVN) distribution (or, equivalently, integrating the MVN density) not only for a range of correlation or covariance structures, but also for a number of dimensions (i.e., variables) that can span several orders of magnitude. In applications for which only one or a few instances of the distribution, and of low dimensionality (\(n\,\hbox to 0.0 pt{\lower 3.0 pt\hbox{$\mathchar 0\sim$}}\raise 1.0 pt\hbox{$<$}\,10\)), must be estimated, conventional numerical methods based on, e.g., Newton-Cotes formulae, Gaussian quadrature and orthogonal polynomials, or tetrachoric series, may offer satisfactory combinations of computational speed and estimation precision.  
Increasingly, however, statistical analysis of large datasets requires many evaluations of very high-dimensional MVN distributions--often as an incidental part of some larger analysis--and places severe demands on the requisite speed and accuracy of numerical methods. We confront the need to estimate the high-dimensional MVN integral in statistical genetics, and particularly in genetic analyses of extended pedigrees (i.e., large, multigenerational collections of related individuals). A typical exercise is variance component analysis of a discrete trait (e.g., a qualitative or categorical measurement of some disease or other condition of interest) under a liability threshold model [1; 2; 3]. Maximum-likelihood estimation of the model parameters in such an application can easily require tens or hundreds of evaluations of the MVN distribution for which \(n\approx\) 100-1000 or greater [4; 5; 6; 7], and situations in which \(n\approx\) 10,000 are not unrealistic.  
In such problems the dimensionality of the model distribution is determined by the product of the total number of individuals in the pedigree(s) to be analyzed and the number of discrete phenotypes jointly analyzed [1; 8]. For univariate traits studied in small pedigrees, such as sibships (sets of individuals born to the same parents) and nuclear families(sibships and their parents), the dimensionality is typically small (\(n\approx 2\)-\(10\)), but analysis of multivariate phenotypes in large extended pedigrees routinely necessitates estimation of mvn distributions for which \(n\) can easily reach several thousand [2; 3; 7]. A single variance component-based linkage analysis of a univariate discrete phenotype in a set of extended pedigrees involves estimating these high-dimensional mvn distributions at hundreds of locations in the genome [3; 9; 10]. In these numerically-intensive applications, estimation of the mvn distribution represents the main computational bottleneck, and the performance of algorithms for estimation of the mvn distribution is of paramount importance.  
Here we report the results of a simulation-based comparison of the performance of two algorithms for estimation of the high-dimensional mvn distribution, the widely-used Mendell-Elston (me) approximation [1; 8; 11; 12] and the Genz Monte Carlo (mc) procedure [13; 14]. Each of these methods is well known, but previous studies have not investigated their properties for very large numbers of dimensions. Using conventional numerical desiderata of estimation accuracy, execution time, and computational efficiency, we examine the performance of these algorithms and identify aspects of the overall mvn context to which each method is particularly well suited.  
In Section 2 we give a brief overview of techniques for estimating the mvn distribution. The me and Genz mc algorithms are reviewed in Section 3. Procedures for exercising the algorithms and comparing their performance are described in Section 4, and results of the comparisons are presented in Section 5. In Section 6 we consider the interpretation and broader implications of our results for future applications.