Skip to content

Double Infinite Integration

4 messages · Aya Anas, S Ellison, David Winsemius +1 more

#
.....
The function z() you provided is not doing anything; it is defined in your post as z<-function(x,y){}

Surprised you didn't get at least NULL, though.








*******************************************************************
This email and any attachments are confidential. Any use...{{dropped:8}}
#
On Dec 4, 2013, at 12:07 PM, Aya Anas wrote:

            
Since f(x) does not depend on y, could you not do:

Int( f(x) * Int( Int( g(y)*z(x,y).dy)).dx)   # not R code

Then use the standard approaches described in several postings to R-help over the years. I'm not a mathematician so anyone is free to correct this.
What we see above is empty space where a body of a function should be.
Surely you didn't enter that!
Perhaps using a more modern tool would be quicker than learning how to do multiple integration with `integrate()`:

http://cran.r-project.org/web/packages/cubature/index.html
Rather that send extraneous exclamation marks, please read the Posting Guide and learn how to get your mail client to send plain text so that the full R code can flow freely onto our devices from yours.
No, no, no.
\/\/\/\/\/\/\/
/\/\/\/\/\/\/\

  
    
#
Aya Anas <aanas <at> feps.edu.eg> writes:
There is a saying: Don't ask Can this be done in R?, ask How is it done?

Extracting function f(x) from the inner integral may not always be the best 
idea. And applying package 'cubature' will not work as adaptIntegrate() does 
not really handle non-finite interval limits.

As an example, let us assume the functions are

    f <- function(x) x
    g <- function(y) y^2
    h <- function(x, y) exp(-(x^2+y^2))

Define a function that calculates the inner integral:

    F1 <- function(x) {
        fun <- function(y) f(x) * g(y) * h(x, y)
        integrate(fun, -Inf, Inf)$value
    }
    F1 <- Vectorize(F1)  # requested when using integrate()

We have to check that integrate() is indeed capable of computing this integrand 
over an infinite interval.

    F1(c(0:4))           # looks good
    ## [1] 0.000000e+00 3.260247e-01 3.246362e-02 3.281077e-04 3.989274e-07

Now integrate this function over the second (infinite) interval.
    
    integrate(F1, 0, Inf)
    ## 0.4431135 with absolute error < 2.4e-06

Correct, as the integral is equal to sqrt(pi)/4 ~ 0.44311346...

If we extract f(x) from the inner integral the value of the integral and the
computation times will be the same, but the overall handling will be slightly 
more complicated.