[R] two dimensional integration
Berend Hasselman
bhh at xs4all.nl
Sun Feb 17 09:41:54 CET 2013
On 16-02-2013, at 18:01, julia cafnik <julia.cafnik at gmail.com> 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]<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
More information about the R-help
mailing list