# 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?

• I don’t know the answer to your question but splinefun can compute derivatives and presumably you could use this as a starting point to fit the data using some inverse methods? I’m interested to learn the solution to this. Sep 27, 2019 at 16:10
• This problem was solved by Maurice Cox in his 1972 paper. I don't know if R supports it, but the search term is "Hermite splines". Sep 27, 2019 at 16:42
• @DavidLeBauer this is what i'm currently doing. I formalised a optimisation problem that fits a number of points such that the spline and it's derivatives approximate the data. But a more direct method would be great.
– dani
Sep 27, 2019 at 17:33
• A quite standard approach is through Kalman filtering. The (unobservable) state contains the exact derivatives and the observations are noisy versions of these. For instance, the model for a cubic spline roughly tells that the second order derivative is a (continuous time) white noise, but a higher-order model can be used as well. You will have to describe the measurement noise depending on the order of derivation for the current observation. Three noise variances (to be estimated) can be enough in a first approach.
– Yves
Oct 1, 2019 at 14:18
• what is measurement error on derivatives? is it much higher than position? also in your plot why arent points aligned? what's x axis? Oct 2, 2019 at 20:47

We will describe how a spline can be used through Kalman Filtering (KF) techniques in relation with a State-Space Model (SSM). The fact that some spline models can be represented by SSM and computed with KF was revealed by C.F. Ansley and R. Kohn in the years 1980-1990. The estimated function and its derivatives are the expectations of the state conditional on the observations. These estimates are computed by using a fixed interval smoothing, a routine task when using a SSM.

For the sake of simplicity, assume that the observations are made at times $$t_1 < t_2 < \dots < t_n$$ and that the observation number $$k$$ at $$t_k$$ involves only one derivative with order $$d_k$$ in $$\{0,\,1,\,2\}$$. The observation part of the model writes as $$\tag{O1} y(t_k) = f^{[d_k]}(t_k) + \varepsilon(t_k)$$ where $$f(t)$$ denotes the unobserved true function and $$\varepsilon(t_k)$$ is a Gaussian error with variance $$H(t_k)$$ depending on the derivation order $$d_k$$. The (continuous time) transition equation takes the general form $$\tag{T1} \frac{\text{d}}{\text{d}t}\boldsymbol{\alpha}(t) = \mathbf{A} \boldsymbol{\alpha}(t) + \boldsymbol{\eta}(t)$$ where $$\boldsymbol{\alpha}(t)$$ is the unobserved state vector and $$\boldsymbol{\eta}(t)$$ is a Gaussian white noise with covariance $$\mathbf{Q}$$, assumed to be independent of the observation noise r.vs $$\varepsilon(t_k)$$. In order to describe a spline, we consider a state obtained by stacking the $$m$$ first derivatives, i.e. $$\boldsymbol{\alpha}(t) := [f(t),\, f^{[1]}(t), \, \dots,\, f^{[m-1]}(t)]^\top$$. The transition is $$\begin{bmatrix} f^{[1]}(t) \\ f^{[2]}(t) \\ \vdots \\ f^{[m-1]}(t) \\ f^{[m]}(t) \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 & &\\ 0 & 0 & 1 & & \\ \vdots & & & \ddots &\\ & & & & 1\\ 0 & \dots & & & 0 \end{bmatrix} \begin{bmatrix} f(t) \\ f^{[1]}(t) \\ \vdots \\ f^{[m-2]}(t)\\ f^{[m-1]}(t) \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ \vdots\\ 0 \\ \eta(t) \end{bmatrix}$$ and we then get a polynomial spline with order $$2m$$ (and degree $$2m-1$$). While $$m=2$$ corresponds to the usual cubic spline, a higher order will be required to use derivatives with order $$>1$$. In order to stick to a classical SSM formalism we can rewrite (O1) as $$\tag{O2} y(t_k) = \mathbf{Z}(t_k) \boldsymbol{\alpha}(t_k) + \varepsilon(t_k),$$ where the observation matrix $$\mathbf{Z}(t_k)$$ picks the suitable derivative in $$\boldsymbol{\alpha}(t_k)$$ and the variance $$H(t_k)$$ of $$\varepsilon(t_k)$$ is chosen depending on $$d_k$$. So $$\mathbf{Z}(t_k) = \mathbf{Z}^\star_{d_k + 1}$$ where $$\mathbf{Z}^\star_1 := [1,\,0,\,\dots,\,0]$$, $$\mathbf{Z}^\star_2 := [0,\,1,\,\dots\,0]$$ and $$\mathbf{Z}^\star_3 := [0,\,0,\,1, 0,\,\dots]$$. Similarly $$H(t_k) = H^\star_{d_k+1}$$ for three variances $$H^\star_1$$, $$H^\star_2$$, and $$H^\star_3$$.

Although the transition is in continuous time, the KF is actually a standard discrete time one. Indeed, we will in practice focus on times $$t$$ where we have an observation, or where we want to estimate the derivatives. We can take the set $$\{t_k\}$$ to be the union of these two sets of times and assume that the observation at $$t_k$$ can be missing: this allows to estimate the $$m$$ derivatives at any time $$t_k$$ regardless of the existence of an observation. There remains to derive the discrete SSM.

We will use indices for discrete times, writing $$\boldsymbol{\alpha}_k$$ for $$\boldsymbol{\alpha}(t_k)$$ and so on. The discrete-time SSM takes the form \begin{align*} \tag{DT} \boldsymbol{\alpha}_{k+1} &= \mathbf{T}_k \,\boldsymbol{\alpha}_{k} + \boldsymbol{\eta}^\star_{k}\\ y_k &= \mathbf{Z}_k\boldsymbol{\alpha}_k + \varepsilon_k \end{align*} where the matrices $$\mathbf{T}_k$$ and $$\mathbf{Q}_k^\star := \text{Var}(\boldsymbol{\eta}_k^\star)$$ are derived from (T1) and (O2) while the variance of $$\varepsilon_k$$ is given by $$H_k=H^\star_{d_k+1}$$ provided that $$y_k$$ is not missing. Using some algebra we can find the transition matrix for the discrete-time SSM $$\mathbf{T}_k = \exp\left\{ \delta_k \mathbf{A} \right\} = \begin{bmatrix} 1 & \frac{\delta_k^1}{1!} & \frac{\delta_k^2}{2!} & \dots & \frac{\delta_k^{m-1}}{(m-1)!}\\ 0 & 1 & \frac{\delta_k^1}{1!} & & \\ \vdots & & & \ddots &\\ & & & & \frac{\delta_k^1}{1!}\\ 0 & \dots & & & 1 \end{bmatrix}, \qquad$$
where $$\delta_k:= t_{k+1} - t_{k}$$ for $$k. Similarly the covariance matrix $$\mathbf{Q}^\star_k = \text{Var} (\boldsymbol{\eta}_k^\star)$$ for the discrete-time SSM can be given as $$\mathbf{Q}^\star_k= \sigma_\eta^2 \, \left[\frac{\delta_k^{2m-i-j+1}}{(m-i)!(m-j)! (2m-i-j+1)}\right]_{i,j}$$ where the indices $$i$$ and $$j$$ are between $$1$$ and $$m$$.

Now to carry over the computation in R we need a package devoted to KF and accepting time-varying models; the CRAN package KFAS seems a good option. We can write R functions to compute the matrices $$\mathbf{T}_k$$ and $$\mathbf{Q}^\star_k$$ from the vector of times $$t_k$$ in order to encode the SSM (DT). In the notations used by the package, a matrix $$\mathbf{R}_k$$ comes to multiply the noise $$\boldsymbol{\eta}^\star_k$$ in the transition equation of (DT): we take it here to be the identity $$\mathbf{I}_m$$. Also note that a diffuse initial covariance must be used here.

EDIT The $$\mathbf{Q}^\star$$ as initially written was wrong. Fixed (aslo in R code and image).

C.F. Ansley and R. Kohn (1986) "On the Equivalence of Two Stochastic Approaches to Spline Smoothing" J. Appl. Probab., 23, pp. 391–405

R. Kohn and C.F. Ansley (1987) "A New Algorithm for Spline Smoothing Based on Smoothing a Stochastic Process" SIAM J. Sci. and Stat. Comput., 8(1), pp. 33–48

J. Helske (2017). "KFAS: Exponential Family State Space Models in R" J. Stat. Soft., 78(10), pp. 1-39

smoothWithDer <- function(t, y, d, m = 3,
Hstar = c(3, 0.2, 0.1)^2, sigma2eta = 1.0^2) {

## define the SSM matrices, depending on 'delta_k' or on 'd_k'
Tfun <- function(delta) {
mat <-  matrix(0, nrow = m, ncol = m)
for (i in 0:(m-1)) {
mat[col(mat) == row(mat) + i] <- delta^i / gamma(i + 1)
}
mat
}
Qfun <- function(delta) {
im <- (m - 1):0
x <- delta^im / gamma(im + 1)
mat <- outer(X = x, Y = x, FUN = "*")
im2 <- outer(im, im, FUN = "+")
sigma2eta * mat * delta / (im2 + 1)
}
Zfun <-  function(d) {
Z <- matrix(0.0, nrow = 1, ncol = m)
Z[1, d + 1] <- 1.0
Z
}
Hfun <- function(d) ifelse(d >= 0, Hstar[d + 1], 0.0)
Rfun <- function() diag(x = 1.0, nrow = m)

## define arrays by stacking the SSM matrices. We need one more
## 'delta' at the end of the series
n <- length(t)
delta <-  diff(t)
delta <- c(delta, mean(delta))

Ta <- Qa <- array(0.0, dim = c(m, m, n))
Za <- array(0.0, dim = c(1, m, n))
Ha <- array(0.0, dim = c(1, 1, n))
Ra <-  array(0.0, dim = c(m, m, n))

for (k in 1:n) {
Ta[ , , k] <- Tfun(delta[k])
Qa[ , , k] <- Qfun(delta[k])
Za[ , , k] <- Zfun(d[k])
Ha[ , , k] <- Hfun(d[k])
Ra[ , , k] <- Rfun()
}

require(KFAS)
## define the SSM and perform Kalman Filtering and smoothing
mod <- SSModel(y ~ SSMcustom(Z = Za, T = Ta, R = Ra, Q = Qa, n = n,
P1 = matrix(0, nrow = m, ncol = m),
P1inf = diag(1.0, nrow = m),
state_names = paste0("d", 0:(m-1))) - 1)
out <- KFS(mod, smoothing = "state")
list(t = t, filtered = out$$att, smoothed = out$$alphahat)

}

## An example function as in OP
f <- function(t, d = rep(0, length = length(t))) {
f <- rep(NA, length(t))
if (any(ind <- (d == 0))) f[ind] <- 2.0 + t[ind] - 0.5 * t[ind]^2
if (any(ind <- (d == 1))) f[ind] <- 1.0 - t[ind]
if (any(ind <- (d == 2))) f[ind] <- -1.0
f
}

set.seed(123)
n <-  100
t <- seq(from = 0, to = 10, length = n)
Hstar <- c(3, 0.4, 0.2)^2
sigma2eta <- 1.0

fTrue <- cbind(d0 = f(t), d1 = f(t, d = 1), d2 = f(t, d = 2))

## ============================================================================
## use a derivative index of -1 to indicate non-observed values, where
## 'y' will be NA
##
## [RUN #0]  no derivative  m = 2 (cubic spline)
## ============================================================================
d0 <- sample(c(-1, 0), size = n, replace = TRUE, prob = c(0.7, 0.3))
ft0 <-  f(t, d0)
## add noise picking the right sd
y0 <- ft0 + rnorm(n = n, sd = c(0.0, sqrt(Hstar))[d0 + 2])
res0 <- smoothWithDer(t, y0, d0, m = 2, Hstar = Hstar)

## ============================================================================
## [RUN #1] Only first order derivative: we can take m = 2 (cubic spline)
## ============================================================================
d1 <- sample(c(-1, 0:1), size = n, replace = TRUE, prob = c(0.7, 0.15, 0.15))
ft1 <-  f(t, d1)
y1 <- ft1 + rnorm(n = n, sd = c(0.0, sqrt(Hstar))[d1 + 2])
res1 <- smoothWithDer(t, y1, d1, m = 2, Hstar = Hstar)

## ============================================================================
## [RUN #2] First and second order derivative: we can take m = 3
## (quintic spline)
## ============================================================================
d2 <- sample(c(-1, 0:2), size = n, replace = TRUE, prob = c(0.7, 0.1, 0.1, 0.1))
ft2 <-  f(t, d2)
y2 <- ft2 + rnorm(n = n, sd = c(0.0, sqrt(Hstar))[d2 + 2])
res2 <- smoothWithDer(t, y2, d2, m = 3, Hstar = Hstar)

## plots : a ggplot with facets would be better here.
for (run in 0:2) {
resrun <- get(paste0("res", run))
drun <- get(paste0("d", run))
yrun <- get(paste0("y", run))
matplot(t, resrun$smoothed, pch = 16, cex = 0.7, ylab = "", xlab = "") matlines(t, fTrue, lwd = 2, lty = 1) for (dv in 0:2) { points(t[drun == dv], yrun[drun == dv], cex = 1.2, pch = 22, lwd = 2, bg = "white", col = dv + 1) } title(main = sprintf("run %d. Dots = smooothed, lines = true, square = obs", run)) legend("bottomleft", col = 1:3, legend = c("d0", "d1", "d2"), lty = 1) } • thank you for your answer. I'm very interested in it. Currently, you don't allow to have the value of f and its derivative to be used at some certain t. Is how is it possible to use all information? Again, merci for your answer. – dani Oct 3, 2019 at 19:48 • My reading is that everything below T1 is about using multiple derivatives in the same inference procedure. Yves can confirm though. Oct 3, 2019 at 22:32 • Indeed, you can use say$o_k >1$derivatives for one$t_k$: the observation$\mathbf{y}_k$is then a vector and$\mathbf{Z}_k$has$o_k$rows picking the wanted derivatives. I am sure that a common$o>1$works with KFAS, but by using NAs it may be possible to have a time varying$o$as well. – Yves Oct 4, 2019 at 5:28 • @Yves Do I understand you correctly: If I have the first and second derivative at point in t_k, then the Z_k looks like this: matrix(c(0,0,0, 0,1,0, 0,0,1), nrow=length(d_k), ncol=m, byrow = T). So, overall it would be a cube of dimension 'highest derivative' * 'spline degree' * '# of time steps' – dani Oct 7, 2019 at 19:31 • Yes @dani, nearly: the number of rows for all the$\mathbf{Z}_k$matrices is$\text{max}_k\{d_k + 1\}$i.e.$3$in the example. This is the highest derivative order plus one. Also, the degree of the spline is$2m-1$, not$m$. In your example since you do not observe the derivative of order$0$(the function itself) it should be set to NA in the observations and you could drop the first row as well. However, I suspect that in this specific case the problem is ill-posed the SSM may not be observable. – Yves Oct 8, 2019 at 5:29 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: 1. The common distribution of the $$\varepsilon_i$$ has a finite variance. 2. 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)

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!