# How to double integrate a function in R

 Classic List Threaded
4 messages
Reply | Threaded
Open this post in threaded view
|

## How to double integrate a function in R

 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

## Re: How to double integrate a function in R

 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]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]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)   { (x0)* (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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

## Re: How to double integrate a function in R

 On Jul 26, 2013, at 9:33 AM, David Winsemius wrote: > fun2 <- function(x,y)   { (x0)* (1/(2*pi))*exp(-y/2)* sqrt((x/(y-x)))} > persp(outer(0:5, 0:6, fun2) ) There does seem to be some potential pathology at the edges of the range, Restricting it to x >= 0.03 removes most of that concern. fun2 <- function(x,y)   { (x0)* (1/(2*pi))*exp(-y/2)* sqrt((x/(y-x)))} persp(outer(seq(0.01,5,by=.01), seq(.02,6,by=.01), fun2) ,ticktype="detailed") > fun <- function(x)   { (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.03,0.03), upper =c(5, 6), tol=1e-2 ) \$integral [1] 0.7605703 \$error [1] 0.00760384 \$functionEvaluations [1] 190859 \$returnCode [1] 0 I tried decreasing the tolerance to 1e-3 but the wait exceeds the patience I have allocated to the problem. -- David Winsemius Alameda, CA, USA ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

## Re: How to double integrate a function in R

 In reply to this post by Tiago Pereira Tiago V. Pereira 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.