[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
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")
curve(f2(x), add = TRUE, col=2)
curve(f3(x), add = TRUE, col=3)
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)'
directly instead of 'splinefun(x$lambda, x$E)'.
Martin Maechler, ETH Zurich
More information about the R-help
mailing list