D1D2: 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.

Value

a list with several components,

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

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

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

spar

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

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

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.

Examples

Run this code

# NOT RUN {
 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")])
 if(is.R()) plot(cos, 0, 10, ylim = c(-1.5,1.5), lwd=2) else { # Splus
   xx <- seq(0,10, len=201); plot(xx, cos(xx), type = 'l', ylim = c(-1.5,1.5))
 }
 title(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)
# }