[R] bivariate normal distribution in likelihood

Peter Dalgaard p.dalgaard at biostat.ku.dk
Tue May 1 11:38:30 CEST 2007

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)).

More information about the R-help mailing list