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. ================================================================ Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail: [EMAIL PROTECTED] 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@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.