[R] bivariate normal distribution in likelihood

[Prof Brian Ripley]

On Tue, 1 May 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).

Which is no worse than other general schemes in high dimensions or 
without smoothness assumptions on the integrand.

> Stratified random sampling would be better, as it converges on the
> order of 1/N.

Your reference for this result, please.  (As stated it is untrue, so I 
presume the stratification scheme depends on N and there are smoothness 
assumptions.)  (BTW, the reference for the results I quote is 'Stochastic 
Simulation' (1987).)

We have not been told 'me' and 'sig', and depending on their values it is 
quite possible that importance sampling would do a great deal better than 
sampling from the specified bivariate normal.

> Even better is product Gauss-Hermite quadrature which will give a
> very accurate answer with a few dozen points.

This is a correlated bivariate normal, and product quadrature methods can 
be arbitrarily bad for such integrals (as people find out for mixed linear 
models).

What you can do is transform to a pair of uncorrelated normals, and for a 
set of the form A as given this transforms to a similar form.  And for 
that an iterated integral can be done easily as the inner integral over y 
will just be a call to pnorm.

Specifically, there is another normal z such that x and z are independent 
and y = z + a*x for some a.  Then A = {x > 0 & z < (3-x)*x + 10} can be 
integrated over z conditional on x and then over x.


> ================================================================
> 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"
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595