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