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

```