[R] integration problem

[Martin Maechler]

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