Here is another approach using `integrate': fvec = function(x, y) sapply(x, function(z, y) exp(-0.549451*(z^2+y^2-0.6*z*y)), y=y) gvec = function(x) sapply(x, function(y) integrate(fvec, lower=-1.51, upper=2.696, subdivisions=10000, rel.tol=1.e-08, y=y)$val)
> integrate(gvec, lower=1.98, upper=3.54, subdivisions=10000, rel.tol=1.e-06) 0.1376102 with absolute error < 1.5e-15 > Ravi. ____________________________________________________________________ Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu ----- Original Message ----- From: Hans W Borchers <hwborch...@googlemail.com> Date: Thursday, July 1, 2010 2:47 pm Subject: Re: [R] Double Integration To: r-h...@stat.math.ethz.ch > Sarah Sanchez <sarah_sanchez09 <at> yahoo.com> writes: > > > > Dear R helpers > > > > I am working on the Bi-variate Normal distribution probabilities. > > I need to double integrate the following function > > (actually simplified form of bivariate normal distribution) > > > > f(x, y) = exp [ - 0.549451 * (x^2 + y^2 - 0.6 * x * y) ] > > > > where 2.696 < x < 3.54 and -1.51 < y < 1.98 > > > > I need to solve something like > > > > INTEGRATE (2.696 to 3.54) dx INTEGRATE [(-1.51 to 1.98)] f(x, y) dy > > > > I have referred to stats::integrate but it deals with only one > variable. > > > > This example appears in Internal Credit Risk Model by Michael Ong > > (page no. 160). > > There has been the same request last month, the answer is still valid: > > library(cubature) > fun <- function(x) exp(-0.549451*(x[1]^2+x[2]^2-0.6*x[1]*x[2])) > adaptIntegrate(fun, c(-1.51,1.98), c(2.696,3.54), tol=1e-8) > > # $integral > # [1] 0.1376102 > # $error > # [1] 1.372980e-09 > # ... > > Hans Werner > > > > > Kindly guide. > > > > Regards > > > > Sarah > > > > ______________________________________________ > R-help@r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org 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.
