R help - May 2007 - bivariate normal distribution in likelihood

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.

Deepankar

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"