Tiago V. Pereira <tiago.pereira <at> mbe.bio.br> writes: > I am trying to double integrate the following expression: > > # expression > (1/(2*pi))*exp(-y2/2)*sqrt((y1/(y2-y1))) > > for y2>y1>0. > > I am trying the following approach > > # first attempt > > library(cubature) > fun <- function(x) { (1/(2*pi))*exp(-x[2]/2)*sqrt((x[1]/(x[2]-x[1])))} > adaptIntegrate(fun, lower = c(0,0), upper =c(5, 6), tol=1e-8) > > However, I don't know how to constrain the integration so that y2>y1>0. > > Any ideas? > Tiago
You could use integral2() in package 'pracma'. It implements the "TwoD" algorithm and has the following properties: (1) The boundaries of the second variable y can be functions of the first variable x; (2) it can handle singularities on the boundaries (to a certain extent). > library(pracma) > fun <- function(y1, y2) (1/(2*pi))*exp(-y2/2)*sqrt((y1/(y2-y1))) > integral2(fun, 0, 5, function(x) x, 6, singular=TRUE) $Q [1] 0.7706771 $error [1] 7.890093e-11 The relative error is a bit optimistic, the absolute error here is < 0.5e-6. The computation time is 0.025 seconds. Hans Werner ______________________________________________ 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.