Re: st: calculation of bivariate normal probabilities

[Arnold Kester <arnold.kester@stat.unimaas.nl>]

Op 05/25/2005 08:18 AM schreef Jenkins S P:
> Stata 9 provides the function -binormal(h,k,r)-, which returns the joint 
> cumulative distribution of the bivariate normal with correlation r; 
> cumulative over (-inf,h]x(-inf,k]. [The equivalent function in version 
> 8.2 was -binorm(h,k,r)-.)
>
> I am not 100% confident about calculation of bivariate normal 
> probabilities once we move outside the standard bivariate normal 
> distribution (or indeed standard univariate normal distribution), and 
> seek advice on this.
>
> My understanding is that -binormal()- evaluates the double integral
>
>         INT(-inf,h] INT(-inf,k) {  f(x1, x2)dx1dx2 }
>
> where f(x1,x2) refers to the standard bivariate normal density with 
> variates X1, X2 each having mean zero, unit variance, and correlation of 
> r. (Put another way, the correlation matrix R has ones on the leading 
> diagonal, and r in the off-diagonal cells)
>
> By contrast, I want to evaluate
>
>     INT(-inf,p] INT(-inf,q) { g(x1,x2)dx1 dx2 }
>
> where g(x1,x2) refers to a general bivariate normal density with 
> variates X1, X2 having means m1 and m2, variances v1, v2, and covariance 
> c (the covariance matrix is C).
>
> To use -binormal()- for the calculation, my inclination was to derive 
> the correlation matrix R from covariance matrix C [-matrix R = corr(C)-] 
> and to extract r from this, and also to set the upper integration points 
> as follows:
>
>     h = (p - m1)/sqrt(v1)
>     q = (q - m2)/sqrt(v2)
>
> Is this correct?
Yes, your transformations standardise the distribution. The covariance 
matrix is standardised by dividing the first row and first column by sd1 
and the second row and column by sd2 (sd=sqrt(v)).

> In fact, I am also seeking an answer for the general multivariate normal 
> case. (If one had a function available to evaluate standard multivariate 
> normal distribution probabilities analogous to -binormal()-, would 
> evaluation for the general case proceed in a similar way? I.e. derive R 
> from C, 'standardize' the upper integration points, and then substitute 
> these into the function.)
The statistical computing program R (cran.r-project.org) has the package 
mvnorm which computes multivariate normal (and T) distribution 
functions, also for higher than two dimensions.

--
Met vriendelijke groet,
Arnold Kester

*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/