Jak dopasować splajn do danych zawierających wartości i 1. / 2. pochodną?

Mam zestaw danych, który zawiera, powiedzmy, kilka pomiarów pozycji, prędkości i przyspieszenia. Wszystkie pochodzą z tego samego „biegu”. Mógłbym zbudować układ liniowy i dopasować wielomian do wszystkich tych pomiarów.

Ale czy mogę zrobić to samo z splajnami? W jaki sposób można to zrobić?

Oto kilka symulowanych danych, które chciałbym dopasować:

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?

Opiszemy, jak można używać splajnu za pomocą technik filtrowania Kalmana (KF) w odniesieniu do modelu State-Space (SSM). Fakt, że niektóre modele splajnu mogą być reprezentowane przez SSM i obliczane za pomocą KF, został ujawniony przez CF Ansley i R. Kohna w latach 1980-1990. Oszacowana funkcja i jej pochodne są oczekiwaniami państwa zależnymi od obserwacji. Oszacowania te są obliczane przy użyciu wygładzania o ustalonych odstępach czasu , rutynowego zadania podczas korzystania z SSM.

Dla uproszczenia załóżmy, że obserwacje są dokonywane w momentach oraz że liczba obserwacji w obejmuje tylko jedną pochodną z rzędem w . Część obserwacyjna modelu zapisuje jako gdzie oznacza nieobserwowaną prawdziwą funkcję, a jest błędem Gaussa z wariancją zależności od rzędu pochodnej . Równanie przejścia (czas ciągły) przyjmuje postać ogólną $t_1 < t_2 < \dots < t_n$$k$$t_k$d k$d_k$$\{0,\,1,\,2\}$

$$\tag{O1} y(t_k) = f^{[d_k]}(t_k) + \varepsilon(t_k)$$ $f(t)$ε ( t k ) H ( t k ) d k d$\varepsilon(t_k)$$H(t_k)$$d_k$ $$\tag{T1} \frac{\text{d}}{\text{d}t}\boldsymbol{\alpha}(t) = \mathbf{A} \boldsymbol{\alpha}(t) + \boldsymbol{\eta}(t)$$ gdzie jest nieobserwowanym wektorem stanu, a to biały szum Gaussa z kowariancją , zakładając, że jest niezależny od hałas obserwacyjny r.vs . Aby opisać splajn, rozważamy stan uzyskany przez ułożenie pierwszych pochodnych, tj. . Przejście jest [f[1](t)f[2](t)⋮f[m-1](t)f[m](t)]=[010$\boldsymbol{\alpha}(t)$$\boldsymbol{\eta}(t)$$\mathbf{Q}$$\varepsilon(t_k)$$m$$\boldsymbol{\alpha}(t) := [f(t),\, f^{[1]}(t), \, \dots,\, f^{[m-1]}(t)]^\top$ $$\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}$$ 2m2m-1m=2>1 rok ( t k ) a następnie otrzymujemy wielomianowy splajn z rzędem (i stopniem ). Podczas gdy odpowiada zwykłem splajnowi sześciennemu,$2m$$2m-1$$m=2$$>1$. Aby trzymać się klasycznego formalizmu SSM, możemy przepisać (O1) jako gdzie macierz obserwacji wybiera odpowiednią pochodną w a wariancja z jest wybierana w zależności od . Więc gdzie , i . Podobnie $$\tag{O2} y(t_k) = \mathbf{Z}(t_k) \boldsymbol{\alpha}(t_k) + \varepsilon(t_k),$$ $\mathbf{Z}(t_k)$$\boldsymbol{\alpha}(t_k)$$H(t_k)$$\varepsilon(t_k)$$d_k$$\mathbf{Z}(t_k) = \mathbf{Z}^\star_{d_k + 1}$$\mathbf{Z}^\star_1 := [1,\,0,\,\dots,\,0]$$\mathbf{Z}^\star_2 := [0,\,1,\,\dots\,0]$$\mathbf{Z}^\star_3 := [0,\,0,\,1, 0,\,\dots]$$H(t_k) = H^\star_{d_k+1}$ H ⋆ 1 H ⋆ 2 H ⋆ 3dla trzech , i . $H^\star_1$$H^\star_2$$H^\star_3$

Chociaż przejście odbywa się w czasie ciągłym, KF jest w rzeczywistości standardowym czasem dyskretnym . Rzeczywiście, w praktyce skupimy się na czasach których mamy obserwację lub w których chcemy oszacować pochodne. Możemy przyjąć zbiór tych dwóch zbiorów czasów i założyć, że może brakować obserwacji w : pozwala to oszacować pochodne dowolnym momencie niezależnie od istnienia obserwacji. Pozostaje jeszcze ustalić dyskretny SSM.$t$$\{t_k\}$$t_k$$m$$t_k$

Użyjemy indeksów dla dyskretnych czasów, pisząc dla i tak dalej. SSM z czasem dyskretnym ma postać gdzie macierze i pochodzą od (T1) i (O2), podczas gdy wariancja jest podana przez pod warunkiem, że$\boldsymbol{\alpha}_k$$\boldsymbol{\alpha}(t_k)$

$$\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*}$$ $\mathbf{T}_k$$\mathbf{Q}_k^\star := \text{Var}(\boldsymbol{\eta}_k^\star)$$\varepsilon_k$$H_k=H^\star_{d_k+1}$$y_k$Tk=exp{δkA}=[ 1 δ 1 knie brakuje. Za pomocą algebry możemy znaleźć macierz przejścia dla dyskretnego czasu SSM gdzie dla . Podobnie macierz kowariancji dla SSM z czasem dyskretnym można podać jako $$\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$$
$\delta_k:= t_{k+1} - t_{k}$$k<n$$\mathbf{Q}^\star_k = \text{Var} (\boldsymbol{\eta}_k^\star)$ $$\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}$$ ij1m gdzie wskaźniki oraz wynoszą od do .$i$$j$$1$$m$

Teraz, aby przenieść obliczenia w R, potrzebujemy pakietu poświęconego KF i akceptującego modele zmieniające się w czasie; pakiet CRAN KFAS wydaje się dobrą opcją. Możemy napisać funkcje R do obliczenia macierzy i z wektora czasów w celu zakodowania SSM (DT). W używanych przez pakiet macierz przychodzi do pomnożenia szumu w równaniu przejścia (DT): bierzemy to tutaj za tożsamość . Należy również pamiętać, że należy tutaj zastosować rozproszoną kowariancję początkową.$\mathbf{T}_k$$\mathbf{Q}^\star_k$$t_k$$\mathbf{R}_k$$\boldsymbol{\eta}^\star_k$$\mathbf{I}_m$

EDIT jak pierwotnie napisany był w błędzie. Naprawiono (również w kodzie R i obrazie).$\mathbf{Q}^\star$

CF Ansley i R. Kohn (1986) „O równoważności dwóch podejść stochastycznych do wygładzania splajnu” J. Appl. Probab , 23, s. 391–405

R. Kohn i CF Ansley (1987) „Nowy algorytm wygładzania splajnu oparty na wygładzaniu procesu stochastycznego” SIAM J. Sci. i Stat. Comput. , 8 (1), s. 33–48

J. Helske (2017). „KFAS: wykładnicze modele przestrzeni stanów rodzinnych w R” J. Stat. Miękki. , 78 (10), s. 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)
}

Możesz zrobić spektakularnie dobrze ze standardową procedurą najmniejszych kwadratów, pod warunkiem, że masz rozsądne pojęcie o względnych rozmiarach błędów losowych popełnianych dla każdej pochodnej. Nie ma ograniczeń co do liczby pomiarów, które wykonujesz dla każdej wartości - możesz nawet jednocześnie mierzyć różne pochodne dla każdej z nich. Jedynym ograniczeniem w stosowaniu zwykłych najmniejszych kwadratów (OLS) jest zwykłe: zakłada się, że pomiary są niezależne.$x$

Podstawową ideę można najlepiej wyrazić poprzez streszczenie problemu. Twój model wykorzystuje zestaw funkcji (jak każda podstawa splajnu) jako podstawa do przewidywania wartości nieznanej funkcji w punktach Oznacza to, że do oszacowania współczynników dla których każda z kombinacji liniowych przybliżeniu zbliżona do Nazwijmy to (wektorową) przestrzenią kombinacji liniowych$p$$f_j:\mathbb{R}\to\mathbb{R},$ $j=1, 2, \ldots, p$$y_i = f(x_i)$$f$$(x_1, x_2, \ldots, x_n).$$\beta_j$$\sum_j \beta_j f_j(x_i)$$y_i.$$\mathbb F.$

Szczególne w tym problemie jest to, że niekoniecznie obserwujesz$y_i.$ Zamiast tego istnieje zdefiniowany zestaw funkcjonałów liniowych powiązany z danymi. Przypomnij sobie, że funkcja jest „funkcją funkcji”: każdy przypisuje liczbę do dowolnej funkcji Model zakłada, że$\mathcal{L}_i$$\mathcal{L}_i$$\mathcal{L}_i[f]$$f\in \mathbb F.$

$$y_i = \mathcal{L}_i [f] + \sigma_i \varepsilon_i\tag{1}$$

gdzie otrzymują funkcjonały, są znanymi współczynnikami skali, a są niezależnymi i identycznie rozmieszczonymi losowymi zmiennymi.$\mathcal{L}_i$$\sigma_i \gt 0$$\varepsilon_i$

Dwa dodatkowe założenia sprawiają, że OLS ma zastosowanie i ma znaczenie statystyczne:

1. Wspólna dystrybucja ma skończoną wariancję.$\varepsilon_i$

2. Każdy jest funkcją liniową . Funkcjonalny jest liniowy, gdy dla dowolnego elementu i odpowiadających mu liczb$\mathcal{L}_i$$\mathcal L$$f_j\in\mathbb{F}$$\alpha_j,$

   $$\mathcal{L}\left[\sum_j \alpha_j f_j\right] = \sum_j \alpha_j \mathcal{L}\left[f_j\right].$$

(2) pozwala na wyraźniejsze wyrażenie modelu jako$(1)$

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

Cały sens tej redukcji polega na tym, że ponieważ określono wszystkie funkcjonały wszystkie podstawowe funkcje a standardowe odchylenia wartości są liczbami - - są to zwykłe „zmienne” lub „cechy” problemu regresji - a są jedynie wagami (względnymi). Tak więc, w optymalnym sensie twierdzenia Gaussa-Markowa, OLS jest świetną procedurą do zastosowania.$\mathcal{L}_i,$$f_j,$$\sigma_i,$$\mathcal{L}_i[f_j]$$\sigma_i$

Funkcjonaliści biorący udział w pytaniu to:

• Oceń w określonym punkcie Tak zwykle robimy. Jest to liniowe, ponieważ z definicji liniowe kombinacje funkcji są oceniane punktowo.$f$$x:$ $\mathcal{L}[f] = f(x).$

• Oceń pochodną w określonym punkcie Jest to liniowe, ponieważ różnicowanie jest liniowe.$f^\prime$$x:$ $\mathcal{L}[f] = f^\prime(x).$

• Oceń drugą pochodną w określonym punkcie$f^{\prime \prime}$$x:$ $\mathcal{L}[f] = f^{\prime \prime}(x).$

---

Ok, jak dobrze działa to podejście? Jak zwykle przestudiujemy resztki porównując dopasowane wartości z wartościami zaobserwowanymi. Ponieważ pozycje, prędkości i przyspieszenia są w różnych jednostkach, należy je narysować na osobnych osiach.$\hat y_i - y_i$$\hat y_i$

Górny rząd używa krzywych do wykreślenia i pierwszych dwóch pochodnych. Odpowiednie punkty danych wykreślono na krzywych: obserwowane wartości po lewej stronie, obserwowane pochodne po środku i obserwowane drugie pochodne po prawej.$\hat y$

Dolny rząd przedstawia odpowiednie wartości resztkowe. Jak zwykle szukamy braku jakiegokolwiek znaczącego związku: mamy nadzieję, że wartości rezydualne (ich współrzędne Y) zmieniają się losowo od lewej do prawej, wykazując niezależność i brak trendów.

Gdy wartości danych wytworzono dokładnie tak, jak w odpowiedzi na pytanie (po ustawieniu losowych liczb do 17 za pomocą na odtwarzalność). Zbadałem pasowania za pomocą spacji B-spline generowanych przez funkcję , również jak w pytaniu, dla stopni od 1 do 6. Ta rycina pokazuje wyniki dla stopnia 2, który jest najniższym stopniem (czyli najprostszym modelem) wykazujący niski AIC i dobre zachowanie resztkowe, a także model wskazany przez ANOVA wszystkich sześciu (zagnieżdżonych) modeli.$n=23$set.seed(17)$\mathbb F$Rbs

Dopasowanie jest

$$\hat y = -27.48993 + 2.54078 f_1 + 2.97679 f_2$$

gdzie i są funkcjami podstawowymi splajnu B utworzonymi przez .$f_1$$f_2$bs

Resztki zachowują się dobrze. Pasowania są dobre. Co więcej, podejście to znalazło właściwy model: dane rzeczywiście zostały wygenerowane z funkcji kwadratowej (stopień 2). Ponadto odchylenia standardowe reszt wynoszą mniej więcej odpowiednie rozmiary: 0,11, 0,20 i 0,61 w porównaniu do 0,1, 0,3 i 0,6 użytych do wygenerowania błędów pierwotnych. To dość niesamowite, biorąc pod uwagę fakt, że krzywe te ekstrapolują obserwacje (które nie wykraczają poza ) i używają tak małego zestawu danych ( ).$x=5$$n=23$

Wreszcie, resztki pasowań dla splajnów wyższego stopnia są jakościowo takie same; wprowadzają jedynie niewielkie ulepszenia kosztem stosowania mniej wiarygodnych modeli. W przypadku wystarczająco wysokich stopni zaczynają się gwałtownie oscylować , na przykład dla małych wartości między obserwowanymi wartościami. Aby zilustrować to (złe) zachowanie, oto dopasowanie stopnia 9:$x$

Wreszcie, oto przykład, w którym dokonano wielu obserwacji różnych funkcjonałów liniowych podstawy. Kod generowania tych obserwacji został zmieniony z kodu w pytaniu na

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

---

RKod do przeprowadzania tych obliczeń jest dość ogólne. W szczególności wykorzystuje numeryczne zróżnicowanie do znalezienia pochodnych, tak aby nie było zależne od rodzaju użytego splajnu. Obsługuje różne wartości obserwacje proporcjonalnie do Automatycznie konstruuje i dopasowuje zestaw modeli w pętli. Funkcjonalności liniowe i odchylenia standardowe są zakodowane na stałe. Istnieją trzy z nich, wybrane zgodnie z wartością zmiennej w zestawie danych.$\sigma_i$$1/\sigma_i^2.$$\mathcal{L}_i$$\sigma_i$type

Jako przykłady tego, w jaki sposób można użyć pasowań, coda drukuje streszczenia, listę ich AIC i ANOVA wszystkich z nich.

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

Przede wszystkim chcę podziękować za postawienie tego pytania. To NAPRAWDĘ interesujące pytanie. Uwielbiam splajny i fajne rzeczy, które możesz z nimi zrobić. To dało mi wymówkę do przeprowadzenia badań. :-)

BLUF: Krótka odpowiedź brzmi „nie”. Nie znam żadnej funkcji w R, która zrobiłaby to dla ciebie automatycznie. Długa odpowiedź jest ... znacznie bardziej skomplikowana. Fakt, że próbki i wartości funkcji nie są próbkowane w tym samym miejscu, utrudnia to. A fakt, że nie masz wartości funkcji blisko prawego końca przedziału, może to uniemożliwić.

Zacznijmy od splajnu sześciennego. Biorąc pod uwagę punkty i odpowiadające im drugie pochodne , przechodzący przez nie splajn sześcienny to:$(x_j, y_j)$$z_j$

$$S_j(x) = Ay_j + By_{j+1} + Cz_j + Dz_{j+1}$$ gdzie Łatwo jest sprawdzić, czy , , i . To gwarantuje, że splajn i jego druga pochodna są ciągłe. Jednak w tym momencie nie mamy ciągłej pierwszej pochodnej. Aby wymusić ciągłość pierwszej pochodnej, potrzebujemy następującego ograniczenia: $$\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}$$ $S_j(x_j) = y_j$$S_j(x_{j+1}) = y_{j+1}$$S''_j(x_j) = z_j$$S''_j(x_{j+1}) = z_{j+1}$ $$\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}$$ W klasycznej konfiguracji splajnu sześciennego zakładasz, że masz punkty i używasz równania (wraz z dwoma dodatkowymi ograniczeniami brzegowymi) do rozwiązania dla . Po znajomości splajn jest w pełni określony i można go użyć do interpolacji w dowolnym dowolnym punkcie. Jako dodatkowy bonus równanie zamienia się w macierz tridiagonalną, którą można rozwiązać w czasie liniowym!$(x_j, y_j)$$\eqref{1}$$z_j$$z_j$$\eqref{1}$

OK, załóżmy teraz, że zamiast znać , znasz . Czy możesz użyć równania aby rozwiązać dla ? Z czystego punktu widzenia algebry wydaje się to wykonalne. Istnieje równań i nieznanych, więc ... dlaczego nie? Ale okazuje się, że nie możesz; matryca będzie pojedyncza. I nie powinno to dziwić. Jak można interpolować wartości funkcji, biorąc pod uwagę TYLKO drugie pochodne? Przynajmniej potrzebujesz wartości początkowej, takiej jak równanie różniczkowe.$y_j$$z_j$$\eqref{1}$$y_j$$N$$N$

Co z twoją sytuacją? Niektóre z twoich punktów mają wartości funkcji, a niektóre z twoich pochodnych. Na razie zignorujmy pierwsze pochodne (to rodzaj bałaganu, z którym trzeba sobie poradzić w przypadku splajnu sześciennego). Formalnie niech będzie zbiorem punktów z wartościami funkcji, a będzie zbiorem punktów z drugimi pochodnymi. Nadal mamy równań z niewiadomymi. tylko, że niektóre niewiadome to a niektóre to . Okazuje się, że otrzymasz rozwiązanie, jeśli 0, 1 lub 2 ORAZ lub$(x_i, y_i), i \in \mathcal{I}$$(x_j, z_j), j \in \mathcal{J}$$N$$N$$y_j$$z_j$$\in \mathcal{I}$$N - 3, N - 2$$N - 1 \in \mathcal{I}$. Innymi słowy, jeden z pierwszych trzech punktów musi być wartością funkcji ORAZ jeden z ostatnich trzech punktów musi być wartością funkcji. Oprócz tego ograniczenia możesz wrzucać dowolną liczbę instrumentów pochodnych.

A co z tymi pierwszymi pochodnymi? Z pewnością możliwe jest uwzględnienie pierwszych pochodnych w splajnie. Ale, jak powiedziałem, robi się znacznie bardziej bałagan. Pierwszą pochodną splajnu podaje: Oczywiście, naprawdę interesuje nas tylko pochodna w węzłach, więc możemy to trochę uprościć, oceniając ją na : można dodawać te ograniczenia macierzy otrzymanej z równania

$$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}$$ $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}$$ $\eqref{1}$a powstały splajn będzie miał określone pierwsze pochodne. Ponadto pomoże to w rozwiązaniu problemu pojedynczej macierzy. Dostaniesz rozwiązanie, jeśli ZNAJDZIESZ wartość funkcji lub pierwszą pochodną w pierwszych trzech i ostatnich trzech punktach.

Więc umieściłem to wszystko w jakimś kodzie i oto zdjęcie, które otrzymałem:

Jak widać, wyniki nie są świetne. Jest tak, ponieważ jest to regularny splajn, który musi honorować WSZYSTKIE dane. Ponieważ dane są stochastyczne, naprawdę musimy użyć splajnu regresji. To temat na inny post. Ale jeśli przejdziesz przez matematykę, skończysz optymalizować kwadratową funkcję celu podlegającą ograniczeniom równości liniowej - i istnieje rozwiązanie w formie zamkniętej!