How can I fit a spline to data that contains values and 1st/2nd derivatives? I have a dataset that contains, let's say, some measurements for position, speed and acceleration. All come from the same "run". I could construct a linear system and fit a polynomial to all of those measurements.
But can I do the same with splines? What is an 'R' way of doing this?
Here is some simulated data I would like to fit:
f <- function(x) 2+x-0.5*x^2+rnorm(length(x), mean=0, sd=0.1)
df <- function(x) 1-x+rnorm(length(x), mean=0, sd=0.3)
ddf <- function(x) -1+rnorm(length(x), mean=0, sd=0.6)

x_f <- runif(5, 0, 5)
x_df <- runif(8, 3, 8)
x_ddf <- runif(10, 4, 9)

data <- data.frame(type=rep('f'), x=x_f, y=f(x_f))
data <- rbind(data, data.frame(type=rep('df'), x=x_df, y=df(x_df)))
data <- rbind(data, data.frame(type=rep('ddf'), x=x_ddf, y=ddf(x_ddf)))

library(ggplot2)
ggplot(data, aes(x, y, color=type)) + geom_point()


library(splines)
m <- lm(data$y ~ bs(data$x, degree=6)) # but I want to fit on f, df, ddf. possible?


 A: You can do spectacularly well with a standard least-squares routine, provided you have a reasonable idea of the relative sizes of the random errors made for each derivative.  There is no restriction on the number of measurements you make for each $x$ value--you can even simultaneously measure different derivatives at each one.  The only limitation in the use of Ordinary Least Squares (OLS) is the usual: you assume the measurements are independent.
The basic idea can be most clearly expressed by abstracting the problem.  Your model uses a set of $p$ functions $f_j:\mathbb{R}\to\mathbb{R},$ $j=1, 2, \ldots, p$ (such as any spline basis) as a basis for predicting the values $y_i = f(x_i)$ of an unknown function $f$ at points $(x_1, x_2, \ldots, x_n).$  This means you seek to estimate coefficients $\beta_j$ for which each of the linear combinations $\sum_j \beta_j f_j(x_i)$ acceptably approximates $y_i.$  Let's call this (vector) space of linear combinations $\mathbb F.$
What is special about this problem is that you don't necessarily observe the $y_i.$  Instead, there is a defined set of linear functionals $\mathcal{L}_i$ associated with the data.  Recall that a functional is a "function of a function:" each $\mathcal{L}_i$  assigns a number $\mathcal{L}_i[f]$ to any function $f\in \mathbb F.$  The model posits that
$$y_i = \mathcal{L}_i [f] + \sigma_i \varepsilon_i\tag{1}$$
where the $\mathcal{L}_i$ are given functionals, the $\sigma_i \gt 0$ are known scale factors, and the $\varepsilon_i$ are independent and identically distributed random variables.
Two additional assumptions make OLS applicable and statistically meaningful:


*

*The common distribution of the $\varepsilon_i$ has a finite variance.

*Every $\mathcal{L}_i$ is a linear functional.  A functional $\mathcal L$ is linear when for any elements $f_j\in\mathbb{F}$ and corresponding numbers $\alpha_j,$ $$\mathcal{L}\left[\sum_j \alpha_j f_j\right] = \sum_j \alpha_j \mathcal{L}\left[f_j\right].$$
(2) permits the model $(1)$ to be expressed more explicitly as

$$y_i = \beta_1 \mathcal{L}_i[f_1] + \cdots + \beta_p \mathcal{L}_i[f_p] + \sigma_i \varepsilon_i.$$

The whole point of this reduction is that because you have stipulated all the functionals $\mathcal{L}_i,$ all the basis functions $f_j,$ and the standard deviations $\sigma_i,$ the values $\mathcal{L}_i[f_j]$ are all numbers--these are just the usual "variables" or "features" of a regression problem--and the $\sigma_i$ are merely (relative) weights.  Thus, in the optimal sense of the Gauss-Markov Theorem, OLS is a great procedure to use.
The functionals involved in the question are the following:


*

*Evaluate $f$ at a specified point $x:$ $\mathcal{L}[f] = f(x).$  This is what we usually do.  This is linear because, by definition, linear combinations of functions are evaluated pointwise.

*Evaluate the derivative $f^\prime$ at a specified point $x:$ $\mathcal{L}[f] = f^\prime(x).$  This is linear because differentiation is linear.

*Evaluate the second derivative $f^{\prime \prime}$ at a specified point $x:$ $\mathcal{L}[f] = f^{\prime \prime}(x).$

Okay, how well does this approach work?  As usual, we will study the residuals $\hat y_i - y_i$ comparing the fitted values $\hat y_i$ to the observed values.  Since positions, velocities, and accelerations are all in different units, they ought to be plotted on separate axes.
 
The top row uses curves to graph $\hat y$ and its first two derivatives.  The relevant data points are plotted over the curves: observed values at the left, observed derivatives in the middle, and observed second derivatives at the right.
The bottom row plots the corresponding residuals.  As usual, we are looking for a lack of any appreciable relationship: we hope the residual values (their y-coordinates) vary randomly from left to right, showing independence and no trends.
The $n=23$ data values were generated exactly as in the question (after setting the random number seed to 17 using set.seed(17) for reproducibility).  I explored fits using the B-spline spaces $\mathbb F$ generated by the R function bs, also as in the question, for degrees 1 through 6.  This figure shows the results for degree 2, which is the lowest degree (that is, simplest model) exhibiting a low AIC and good residual behavior, as well as the model indicated by an ANOVA of all six (nested) models.
The fit is
$$\hat y = -27.48993 + 2.54078 f_1 + 2.97679 f_2$$
where $f_1$ and $f_2$ are the B-spline basis functions created by bs.
The residuals behave well.  The fits are good. Moreover, this approach found the correct model: the data indeed were generated from a quadratic function (degree 2).  Furthermore, the standard deviations of the residuals are about the right sizes: 0.11, 0.20, and 0.61 compared to 0.1, 0.3, and 0.6 used to generate the original errors.  That's pretty amazing given that these curves obviously extrapolate the observations (which do not go beyond $x=5$) and use such a small dataset ($n=23$).
Finally, residuals to the fits for higher-degree splines are qualitatively the same; they make only slight improvements at a cost of using less-plausible models.  For sufficiently high degrees, they begin to oscillate wildly for small values of $x$ between the observed values, for instance.  To illustrate this (bad) behavior, here's the degree-9 fit:

Finally, here is an example where multiple observations of various linear functionals of the basis were made.  The code for generating these observations was changed from that in the question to
mult <- 2
x_f <- rep(runif(5, 0, 5), mult)       # Two observations per point
x_df <- rep(runif(8, 3, 8), mult)      # Two derivatives per point
x_ddf <- c(x_df, rep(runif(10, 4, 9))  # Derivative and acceleration per point



The R code for carrying these calculations is rather general.  In particular, it uses numerical differentiation to find the derivatives so that it is not dependent on the type of spline used.  It handles the differing values of $\sigma_i$ by weighting the observations proportionally to $1/\sigma_i^2.$  It automatically constructs and fits a set of models in a loop.  The linear functionals $\mathcal{L}_i$ and the standard deviations $\sigma_i$ are hard-coded.  There are three of each, selected according to the value of the type variable in the dataset.
As examples of how you can use the fits, the coda prints summaries, a list of their AICs, and an ANOVA of them all.
#
# Estimate spline derivatives at points of `x`.
#
d <- function(x, s, order=1) {
  h <- diff(range(x, na.rm=TRUE))
  dh <- h * 1e-4
  lags <- seq(-order, order, length.out=order+1) * dh/2
  b <- choose(order, 0:order) * (-1)^(order:0)
  y <- b %*% matrix(predict(s, c(outer(lags, x, `+`))), nrow=length(lags))
  y <- matrix(y / (dh^order), nrow=length(x))
}
#
# Fit and plot models by degree.
#
data$order <- c(f=0, df=1, ddf=2)[data$type]
k <- max(data$order)
x <- data$x
w <- (c(0.1, 0.3, 0.6)^(-2))[data$order+1] # As specified in the question

fits <- lapply(1:6, function(deg) {
  #
  # Construct a model matrix.
  #
  s <- bs(x, degree=deg, intercept=TRUE)
  X.l <- lapply(seq.int(k+1)-1, function(i) {
    X <- subset(data, order==i)
    Y <- as.data.frame(d(X$x, s, order=i))
    cbind(X, Y)
  })
  X <- do.call("rbind", X.l)
  #
  # Fit WLS models.
  #
  f <- as.formula(paste("y ~ -1 +", paste0("V", 0:deg+1, collapse="+")))
  fit <- lm(f, X, weights=w)
  msr <- tapply(residuals(fit), data$order, function(r) {
    k <- length(r) - 1 - deg
    ifelse(k >= 1, sum(r^2) / k, 1)
  })
  #
  # Compute predicted values along the graphs.
  #
  X.new <- data.frame(x = seq(min(X$x), max(X$x), length.out=101))
  X.new$y.hat <- predict(s, X.new$x) %*% coefficients(fit)
  X.new$Dy.hat <- d(X.new$x, s, 1) %*% coefficients(fit)
  X.new$DDy.hat <- d(X.new$x, s, 2) %*% coefficients(fit)
  X$Residual <- residuals(fit)
  #
  # Return the model.
  #
  fit$msr <- msr
  fit
})
lapply(fits, function(f) sqrt(f$msr))
lapply(fits, summary)
lapply(fits, AIC)
do.call("anova", fits)

A: First of all, I want to thank you for posing this question.  It's a REALLY interesting question.  I love splines and the cool things you can do with them.  And this gave me an excuse to do some research.  :-)
BLUF: The short answer is no.  I don't know of any functionality in R that will do this for you automatically.  The long answer is... much more complicated.  The fact that the derivatives and function values aren't sampled at the same place makes this more difficult.  And the fact that you don't have a function value near the right end of the interval might make it impossible.
Let's start with the cubic spline.  Given points $(x_j, y_j)$ and the corresponding second derivatives $z_j$, the cubic spline passing through them is:
$$
S_j(x) = Ay_j + By_{j+1} + Cz_j + Dz_{j+1}
$$
where
$$
\begin{array}{}
h_j & = & x_{j+1} - x_j \\
A & = & \frac{x_{j+1} - x}{h_j} \\
B & = & 1 - A \\
C & = & \frac{1}{6}(A^3 - A)h_j ^2 \\
D & = & \frac{1}{6}(B^3 - B)h_j ^2
\end{array}
$$
It's pretty straightforward to verify that $S_j(x_j) = y_j$, $S_j(x_{j+1}) = y_{j+1}$, $S''_j(x_j) = z_j$ and $S''_j(x_{j+1}) = z_{j+1}$.  This guarantees that the spline and its second derivative are continuous.  However, at this point, we don't have a continuous first derivative. In order to force the first derivative to be continuous, we need the following constraint:
$$
\frac{6}{h_{j-1}}y_{j-1} - \left( \frac{6}{h_{j-1}} + \frac{6}{h_j} \right) y_j + \frac{6}{h_j}y_{j+1} = h_{j-1} z_{j-1} + 2(h_{j-1} + h_j) z_j + h_j z_{j + 1}  \tag{1}\label{1}
$$
In the classic cubic spline setup, you assume you have the points $(x_j, y_j)$ and use equation \eqref{1} (along with two additional boundary constraints) to solve for the $z_j$.  Once you know the $z_j$, the spline is fully specified and you can use it to interpolate at any arbitrary point.  As an added bonus, equation \eqref{1} turns in to a tridiagonal matrix which can be solved in linear time!
OK, now suppose that, instead of knowing the $y_j$, you know the $z_j$.  Can you use equation \eqref{1} to solve for the $y_j$?  From a pure algebra standpoint, it seems feasible.  There are $N$ equations and $N$ unknowns, so... why not?  But it turns out you can't; the matrix will be singular.  And that should come as no surprise.  How could you possibly interpolate the function values given JUST the second derivatives?  At the very least, you would need an initial value, just like a differential equation.
What about your situation?  Some of your points have function values and some of your points have derivatives.  For the time being, let's ignore the first derivatives (they're kind of a mess to deal with in the cubic spline basis). 
 Formally, let $(x_i, y_i), i \in \mathcal{I}$ be the set of points with function values and $(x_j, z_j), j \in \mathcal{J}$ be the set of points with second derivatives.  We still have $N$ equations with $N$ unknowns.  It's just that some of the unknowns are $y_j$ and some are $z_j$.  It turns out that you will get a solution if 0, 1 or 2 $\in \mathcal{I}$ AND $N - 3, N - 2$ or $N - 1 \in \mathcal{I}$.  In other words, one of the first three points has to be a function values AND one of the last three points has to be a function value.  Other than that constraint, you're free to throw in as many derivatives as you want.
How about those first derivatives?  It's certainly possible to include first derivatives in your spline.  But, like I said, it gets a lot messier.  The first derivative of the spline is given by:
$$
S'_j(x) = \frac{y_{j+1} - y_j}{h_j} - \frac{3A^2 - 1}{6} h_j z_j + \frac{3B^2 - 1}{6} h_j z_{j+1}
$$
Of course, we're only really interested in the derivative at the knots, so we can simplify this a little bit by evaluating it at $x_j$:
$$
S'_j(x_j) = \frac{y_{j+1} - y_j}{h_j} - \frac{1}{3} h_j z_j - \frac{1}{6} h_j z_{j+1}
$$
You can add these constraints to the matrix you get from equation \eqref{1} and the resulting spline will have the specified first derivatives.  In addition, this will help with the singular matrix problem.  You'll get a solution if you have EITHER a function value or a first derivative in the first three and last three points.
So I put that all together in some code and here's the picture I got:

As you can see, the results aren't great.  That's because this is a regular spline which must honor ALL of the data.  Since the data is stochastic, we really need to use a regression spline.  That's a topic for another post.  But if you work through the math, you'll end up optimizing a quadratic objective function subject to linear equality constraints - and there's a closed form solution!
