On 16-02-2013, at 18:01, julia cafnik <[email protected]> wrote:
> Dear R-users, > > I'm wondering how to calculate this double integral in R: > > int_a^b int_c^y g(x, y) dx dy > > where g(x,y) = exp(- alpha (y - x)) * b > A very similar question was asked about nine years ago: http://tolstoy.newcastle.edu.au/R/help/04/10/5831.html The final message in the thread gave a workable answer. In your case define function g (leave out the multiply by b since you can always do that outside the integral). g <- function(x,y,alpha=1) exp(-alpha*(y-x)) Assume a=0, b=1 and c=0. Then following the final message in the thread mentioned above there are two ways of getting an answer: if function g is fully vectorized: integrate( function(y) { sapply(y, function(y) { integrate(function(x) g(x,y),0,y)$value }) },0,1) and if g is not vectorized and only takes scalars as input integrate(function(y) { sapply(y, function(y) { integrate(function(x) { sapply(x, function(x) g(x,y)) },0,y)$value }) },0,1) Both of the approaches result in 0.3678794 with absolute error < 4.1e-15 which corresponds to the exact analytical answer 1/exp(1) (if I didn't make a mistake) You can also do this with package cubature with Gabor's approach (from the mentioned thread) like this library(cubature) h <- function(z) g(z[1],z[2])*(z[1]<z[2]) adaptIntegrate(h,c(0,0),c(1,1),tol=1e-4) but it's also very slow with higher tolerances. > adaptIntegrate(h,c(0,0),c(1,1),tol=1e-4) $integral [1] 0.3678723 $error [1] 3.677682e-05 $functionEvaluations [1] 268413 $returnCode [1] 0 Berend ______________________________________________ [email protected] 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.
