two dimensional integration

Dear R-users,

I'm wondering how to calculate this double integral in R:

int_a^b  int_c^y g(x, y) dx dy

where g(x,y) = exp(- alpha (y - x)) * b

Thanks for answering!

One way would be to install the cubature package.

David Winsemius
Alameda, CA, USA

Hello,

Though I think you should compute that integral symbolically by hand and 
then define a function with the result, maybe package pracma can do what 
you want.
[functions dblquad(9 and quad2d()]


Hope this helps,

Rui Barradas

Em 16-02-2013 17:01, julia cafnik escreveu:

Dear R-users,

I'm wondering how to calculate this double integral in R:

int_a^b  int_c^y g(x, y) dx dy

where g(x,y) = exp(- alpha (y - x)) * b

Thanks for answering!

Cheers,

Alui

	[[alternative HTML version deleted]]

______________________________________________
R-help at r-project.org 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.

Dear R-users,

I'm wondering how to calculate this double integral in R:

int_a^b  int_c^y g(x, y) dx dy

where g(x,y) = exp(- alpha (y - x)) * b

A very similar question was asked about nine years ago: http://tolstoy.newcastle.edu.au/R/help/04/10/5831.html
The final message in the thread gave a workable answer.

In your case define function g (leave out the multiply by b since you can always do that outside the integral).

g <- function(x,y,alpha=1) exp(-alpha*(y-x))

Assume a=0, b=1 and c=0.

Then following the final message in the thread mentioned above there are two ways of getting an answer:

if function g is fully vectorized:

integrate( function(y) {
    sapply(y, function(y) { 
        integrate(function(x) g(x,y),0,y)$value
    })
},0,1)

and if g is not vectorized and only takes scalars as input

integrate(function(y) {
  sapply(y, function(y) {
    integrate(function(x) {
      sapply(x, function(x) g(x,y)) },0,y)$value
  })
},0,1)

Both of the approaches result in

0.3678794 with absolute error < 4.1e-15

which corresponds to the exact analytical answer 1/exp(1) (if I didn't make a mistake)

You can also do this with package cubature with Gabor's approach (from the mentioned thread) like this

library(cubature)              
h <- function(z) g(z[1],z[2])*(z[1]<z[2])
adaptIntegrate(h,c(0,0),c(1,1),tol=1e-4) 

but it's also very slow with higher tolerances.

adaptIntegrate(h,c(0,0),c(1,1),tol=1e-4)

$integral
[1] 0.3678723
$error
[1] 3.677682e-05
$functionEvaluations
[1] 268413
$returnCode
[1] 0


Berend

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20130217/372b0016/attachment.pl>

thank for your help. already solved it.

Show us how.
So that others looking for answers to similar problems in future can find an answer.

Berend

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20130217/144e986e/attachment.pl>