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



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Oct 27 18:19:05 2004

This archive was generated by hypermail 2.1.8 : Wed 03 Nov 2004 - 22:58:25 EST