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