Hello, R users! 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
On Jul 26, 2013, at 8:44 AM, Tiago V. Pereira wrote:> 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.Generally incorporating boundaries is accomplished by multiplying the integrand with logical vectors that encapsulate what are effectively two conditions: Perhaps: fun <- function(x) { (x[1]<x[2])*(x[1]>0)* (1/(2*pi))*exp(-x[2]/2)* sqrt((x[1]/(x[2]-x[1])))} That was taking quite a long time and I interrupted it. There were quite a few warnings of the sort 1: In sqrt((x[1]/(x[2] - x[1]))) : NaNs produced 2: In sqrt((x[1]/(x[2] - x[1]))) : NaNs produced Thinking the NaNs might sabotage the integration process, I added a conditional to the section of that expression that was generating the NaNs. I don't really know whether NaN's are excluded from the summation process in adaptIntegrate: fun <- function(x) { (x[1]<x[2])*(x[1]>0)* (1/(2*pi))*exp(-x[2]/2)* if(x[1]>x[2]){ 0 }else{ sqrt((x[1]/(x[2]-x[1])) )} } adaptIntegrate(fun, lower = c(0,0), upper =c(5, 6) ) I still didn't have the patience to wait for an answer, but I did plot the function: fun2 <- function(x,y) { (x<y)*(x>0)* (1/(2*pi))*exp(-y/2)* sqrt((x/(y-x)))} persp(outer(0:5, 0:6, fun2) ) So at least the function is finite over most of its domain. -- David Winsemius Alameda, CA, USA
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? > TiagoYou 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
I would like to express my gratitude for the great help given by David and Hans regarding my last query. Thank you very much for your time, people. All the best, Tiago --- Hello, R users! 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 -- Tiago V. Pereira, MSc, PhD Center for Studies of the Human Genome Department of Genetics and Evolutionary Biology University of S?o Paulo Rua do Mat?o, 277 CEP 05508-900 S?o Paulo - SP, Brazil