How to double integrate a function in R

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

How to double integrate a function in R

Tiago Pereira
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: How to double integrate a function in R

David Winsemius

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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: How to double integrate a function in R

David Winsemius

On Jul 26, 2013, at 9:33 AM, David Winsemius wrote:

> 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) )

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)   { (x<y)*(x>0)* (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]<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.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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|

Re: How to double integrate a function in R

Hans W Borchers
In reply to this post by Tiago Pereira
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

______________________________________________
[hidden email] 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.