D1D2 {sfsmisc}

Numerical Derivatives of (x,y) Data via Smoothing Splines

Description

Compute numerical derivatives of f() given observations (x,y), using cubic smoothing splines with GCV, see smooth.spline. In other words, estimate f'() and/or f''() for the model

Y_i = f(x_i) + E_i, \ \ i = 1,\dots n,

Usage

D1D2(x, y, xout = x, spar.offset = 0.1384, deriv = 1:2, spl.spar = NULL)

Arguments

x, y

numeric vectors of same length, supposedly from a model y ~ f(x).

xout

abscissa values at which to evaluate the derivatives.

spar.offset

numeric fudge added to the smoothing parameter, see spl.par below.

deriv

integer in 1:2 indicating which derivatives are to be computed.

spl.spar

direct smoothing parameter for smooth.spline. If it is NULL (as per default), the smoothing parameter used will be spar.offset + sp$spar, where sp$spar is the GCV estimated smoothing parameter, see smooth.spline.

Details

It is well known that for derivative estimation, the optimal smoothing parameter is larger (more smoothing) than for the function itself. spar.offset is really just a fudge offset added to the smoothing parameter. Note that in R's implementation of smooth.spline, spar is really on the \log\lambda scale.

When deriv = 1:2 (as per default), both derivatives are estimated with the same smoothing parameter which is suboptimal for the single functions individually. Another possibility is to call D1D2(*, deriv = k) twice with k = 1 and k = 2 and use a larger smoothing parameter for the second derivative.

Value

a list with several components,

x

the abscissae values at which the derivative(s) are evaluated.

D1

if deriv contains 1, estimated values of f'(x_i) where x_i are the values from xout.

D2

if deriv contains 2, estimated values of f''(x_i).

spar

the smoothing parameter used in the (final) smooth.spline call.

df

the equivalent degrees of freedom in that smooth.spline call.

Author(s)

Martin Maechler, in 1992 (for S).

See Also

D2ss which calls smooth.spline twice, first on y, then on the f'(x_i) values; smooth.spline on which it relies completely.

Examples

 set.seed(8840)
 x <- runif(100, 0,10)
 y <- sin(x) + rnorm(100)/4

 op <- par(mfrow = c(2,1))
 plot(x,y)
 lines(ss <- smooth.spline(x,y), col = 4)
 str(ss[c("df", "spar")])
 plot(cos, 0, 10, ylim = c(-1.5,1.5), lwd=2,
      main = expression("Estimating f'() : "~~ frac(d,dx) * sin(x) == cos(x)))
 offs <- c(-0.1, 0, 0.1, 0.2, 0.3)
 i <- 1
 for(off in offs) {
   d12 <- D1D2(x,y, spar.offset = off)
   lines(d12$x, d12$D1, col = i <- i+1)
 }
 legend(2,1.6, c("true cos()",paste("sp.off. = ", format(offs))), lwd=1,
        col = 1:(1+length(offs)), cex = 0.8, bg = NA)
 par(op)