# Re: [R] integration problem

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Wed 27 Oct 2004 - 18:02:34 EST

>>>>> "Sundar" == Sundar Dorai-Raj <sundar.dorai-raj@pdf.com> >>>>> on Tue, 26 Oct 2004 11:19:36 -0500 writes:

Sundar> Christoph Scherber wrote:
>> Dear R users,
>>
>> I have spectral data (say, wavelength vs. extinction coefficient) for
>> which IŽd like to calculate an integral (i.e. the area underneath the
>> curve).
>>
>> Suppose the (artificial) dataset is
>>
>> lambda E
>> 1 2
>> 2 4
>> 3 5
>> 4 8
>> 5 1
>> 6 5
>> 7 4
>> 8 9
>> 9 8
>> 10 2
>>
>>
>>
>> How can I calculate an integral for these values using R?
>>
>> Many thanks for any help!
>> Regards
>>
>> Christoph
>>

```    Sundar> Hi Christoph,
Sundar> You can write some code to do trapezoidal integration or use ?approx
Sundar> in the following manner:

```

Sundar> f <- function(xnew, x, y) approx(x, y, xnew)\$y     Sundar> integrate(f, min(x\$lambda), max(x\$lambda), x = x\$lambda, y = x\$E)

```    Sundar> where `x' is your data above. Using approx is
Sundar> perhaps overkill but it works. A better solution
Sundar> would be to use trapezoids or perhaps piecewise
Sundar> linear integration. I don't know of a package that
Sundar> has the latter approaches off the top of my head but
Sundar> that doesn't mean they doesn't exist somewhere.

```

And then, it mighgt be slightly more satisfactory to use spline() instead of approx() {smooth interpolation instead of linear}. However, this is exactly an example where  approxfun() is much nicer than approx() and  splinefun() is much nicer than spline() :

I here give a script {output in comments}:

x <- data.frame(lambda = 1:10, E = c(2,4,5,8,1, 5,4,9,8,2))

## Sundar's proposal
f <- function(xnew, x, y) approx(x, y, xnew)\$y  integrate(f, min(x\$lambda), max(x\$lambda), x = x\$lambda, y = x\$E)  ##--> 46 with absolute error < 0.0047

## Using the *fun() version

f2 <- approxfun(x\$lambda, x\$E)
## or even
f2 <- approxfun(x) # does work: 2-column data frame or matrix

integrate(f2, min(x\$lambda), max(x\$lambda))   ## same answer as above

f3 <- splinefun(x\$lambda, x\$E)
integrate(f3, min(x\$lambda), max(x\$lambda))   ##--> 46.98437 with absolute error < 0.0046

## Note that you can do more than just integrate them:

plot(x, type = "b")
f4 <- splinefun(x\$lambda, x\$E, method ="natural") curve(f4(x), add = TRUE, col=4)# not much difference

###-----

Further note, that only in R-patched (or R-devel), we can also use 'splinefun(x)'