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-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.