[R] bivariate normal distribution in likelihood

[DEEPANKAR BASU]

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