R help - Jul 2013 - 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

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

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