Numerical Derivatives of (x,y) Data via Smoothing Splines
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,… n,
D1D2(x, y, xout = x, spar.offset = 0.1384, deriv = 1:2, spl.spar = NULL)
x,y |
numeric vectors of same length, supposedly from a model
|
xout |
abscissa values at which to evaluate the derivatives. |
spar.offset |
numeric fudge added to the smoothing parameter,
see |
deriv |
integer in |
spl.spar |
direct smoothing parameter for |
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λ 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.
a list with several components,
x |
the abscissae values at which the derivative(s) are evaluated. |
D1 |
if |
D2 |
if |
spar |
the smoothing parameter used in the (final)
|
df |
the equivalent degrees of freedom in that
|
Martin Maechler, in 1992 (for S).
D2ss which calls smooth.spline twice,
first on y, then on the f'(x_i) values;
smooth.spline on which it relies completely.
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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.