[R] bivariate normal distribution in likelihood

[Peter Dalgaard]

Robert A LaBudde wrote:
> 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.
>   
Or, use the mvtnorm package. It has pmvnorm, which does the integrals 
over rectangular regions. You'll need a pretransformation to use it for 
the problem at hand, though (from (x,y) to (x,y-3x)).