1\name{D2ss}
2\alias{D2ss}
3\alias{D1ss}
4\alias{D1tr}
5\title{Numerical Derivatives of (x,y) Data  (via Smoothing Splines)}
6\description{
7  Compute the numerical first or 2nd derivatives of \eqn{f()} given
8  observations \code{(x[i], y ~= f(x[i]))}.
9
10  \code{D1tr} is the \emph{\bold{tr}ivial} discrete first derivative
11  using simple difference ratios, whereas \code{D1ss} and \code{D2ss}
12  use cubic smoothing splines (see \code{\link[stats]{smooth.spline}})
13  to estimate first or second derivatives, respectively.
14
15  \code{D2ss} first uses \code{smooth.spline} for the first derivative
16  \eqn{f'()} and then applies the same to the predicted values
17  \eqn{\hat f'(t_i)}{f'^(t[i])} (where \eqn{t_i}{t[i]} are the values of
18  \code{xout}) to find   \eqn{\hat f''(t_i)}{f''^(t[i])}.
19}
20\usage{
21D1tr(y, x = 1)
22
23D1ss(x, y, xout = x, spar.offset = 0.1384, spl.spar=NULL)
24D2ss(x, y, xout = x, spar.offset = 0.1384, spl.spar=NULL)
25}
26\arguments{
27  \item{x,y}{numeric vectors of same length, supposedly from a model
28    \code{y ~ f(x)}.  For \code{D1tr()}, \code{x} can have length one
29    and then gets the meaning of \eqn{h = \Delta x}.}
30  \item{xout}{abscissa values at which to evaluate the derivatives.}
31  \item{spar.offset}{numeric fudge added to the smoothing parameter(s),
32    see \code{spl.par} below.  Note that the current default is there
33    for historical reasons only, and we often  would recommend to use
34    \code{spar.offset = 0} instead.}
35  \item{spl.spar}{direct smoothing parameter(s) for \code{smooth.spline}.
36    If it is \code{NULL} (as per default), the smoothing parameter used
37    will be \code{spar.offset + sp$spar}, where \code{sp$spar} is the GCV
38  estimated smoothing parameter for \emph{both} smooths, see
39  \code{\link{smooth.spline}}.}
40}
41\details{
42  It is well known that for derivative estimation, the optimal smoothing
43  parameter is larger (more smoothing needed) than for the function itself.
44  \code{spar.offset} is really just a \emph{fudge} offset added to the
45  smoothing parameters. Note that in \R's implementation of
46  \code{\link{smooth.spline}}, \code{spar} is really on the
47  \eqn{\log\lambda} scale.
48%
49%   When \code{deriv = 1:2} (as per default), both derivatives are
50%   estimated with the \emph{same} smoothing parameter which is suboptimal
51%   for the single functions individually.  Another possibility is to call
52%   \code{D1D2(*, deriv = k)} twice with \code{k = 1} and \code{k = 2} and
53%   use a \emph{larger} smoothing parameter for the second derivative.
54}
55\value{
56  \code{D1tr()} and \code{D1ss()} return a numeric vector of the length
57  of \code{y} or \code{xout}, respectively.
58
59  \code{D2ss()} returns a list with components
60  \item{x}{the abscissae values (= \code{xout}) at which the
61    derivative(s) are evaluated.}
62  \item{y}{estimated values of \eqn{f''(x_i)}.}
63  \item{spl.spar}{numeric vector of length 2, contain the \code{spar}
64    arguments to the two \code{smooth.spline} calls.}
65  \item{spar.offset}{as specified on input (maybe rep()eated to length 2).}
66}
67\author{Martin Maechler, in 1992 (for S).}
68\seealso{\code{\link{D1D2}} which directly uses the 2nd derivative of
69  the smoothing spline;  \code{\link{smooth.spline}}.
70}
71\examples{
72
73## First Derivative  --- spar.off = 0  ok "asymptotically" (?)
74set.seed(330)
75mult.fig(12)
76for(i in 1:12) {
77  x <- runif(500, 0,10); y <- sin(x) + rnorm(500)/4
78  f1 <- D1ss(x=x,y=y, spar.off=0.0)
79  plot(x,f1, ylim = range(c(-1,1,f1)))
80  curve(cos(x), col=3, add= TRUE)
81}
82
83 set.seed(8840)
84 x <- runif(100, 0,10)
85 y <- sin(x) + rnorm(100)/4
86
87 op <- par(mfrow = c(2,1))
88 plot(x,y)
89 lines(ss <- smooth.spline(x,y), col = 4)
90 str(ss[c("df", "spar")])
91 xx <- seq(0,10, len=201)
92 plot(xx, -sin(xx), type = 'l', ylim = c(-1.5,1.5))
93 title(expression("Estimating f''() :  " * frac(d^2,dx^2) * sin(x) == -sin(x)))
94 offs <- c(0.05, 0.1, 0.1348, 0.2)
95 i <- 1
96 for(off in offs) {
97   d12 <- D2ss(x,y, spar.offset = off)
98   lines(d12, col = i <- i+1)
99 }
100 legend(2,1.6, c("true :  -sin(x)",paste("sp.off. = ", format(offs))), lwd=1,
101        col = 1:(1+length(offs)), cex = 0.8, bg = NA)
102 par(op)
103}
104\keyword{smooth}
105