[R] integration problem

Martin Maechler maechler at stat.math.ethz.ch
Wed Oct 27 10:02:34 CEST 2004

```>>>>> "Sundar" == Sundar Dorai-Raj <sundar.dorai-raj at 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
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))

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)'