[R] bivariate normal distribution in likelihood

[Robert A LaBudde]

At 11:32 PM 4/30/2007, Deepankar wrote:
>Hi all,
>
>I am trying to do a maximum likelihood estimation. In my likelihood 
>function, I have to evaluate a double integral of the bivariate 
>normal density over a subset of the points of the plane. For 
>instance, if I denote by y the y co-ordinate and by x, the x 
>co-ordinate then the area over which I have to integrate is the 
>following set of points, A:
>A = {x>=0 & y < 3x+10}
>
>I have used the following code to calculate this double integral:
>
>x <- rmvnorm(100000, mean=me, sigma=sig)
>c <- nrow(x)
>x1 <- x[ ,1]
>x2 <- x[ ,2]
>e1 <- as.numeric(x2 < 3*x1 + 10 & x1>0)
>p1 <- sum(e1)/c
>
>In this code, I provide the mean and covariance while drawing the 
>bivariate random normal variables and get "p1" as the required 
>answer. The problem is that I have to draw at least 100,000 
>bivariate random normals to get a reasonable answer; even then it is 
>not very accurate.
>
>Is there some other way to do the same thing more accurately and 
>more efficiently? For instance, can this be done using the bivariate 
>normal distribution function "pmvnorm"? Also feel free to point our 
>errors if you see one.

Simple random sampling is a poor way to evaluate an integral 
(expectation). It converges on the order of 1/sqrt(N).

Stratified random sampling would be better, as it converges on the 
order of 1/N.

Even better is product Gauss-Hermite quadrature which will give a 
very accurate answer with a few dozen points.

================================================================
Robert A. LaBudde, PhD, PAS, Dpl. ACAFS  e-mail: ral at lcfltd.com
Least Cost Formulations, Ltd.            URL: http://lcfltd.com/
824 Timberlake Drive                     Tel: 757-467-0954
Virginia Beach, VA 23464-3239            Fax: 757-467-2947

"Vere scire est per causas scire"