thank for your help. already solved it.
Cheers,
J.
On Sun, Feb 17, 2013 at 9:41 AM, Berend Hasselman wrote:
>
> On 16-02-2013, at 18:01, julia cafnik wrote:
>
> > 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] 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
>
>
[[alternative HTML version deleted]]