Okresowe splajny w celu dopasowania danych okresowych

W komentarzu do tego pytania użytkownik @whuber wskazał na możliwość używania okresowej wersji splajnów w celu dopasowania danych okresowych. Chciałbym dowiedzieć się więcej o tej metodzie, w szczególności o równaniach definiujących splajny oraz o tym, jak je zaimplementować w praktyce (jestem głównie Rużytkownikiem, ale jeśli zajdzie taka potrzeba, poradzę sobie z MATLAB-em lub Pythonem). Ponadto, ale jest to „miło mieć”, dobrze byłoby wiedzieć o możliwych zaletach / wadach w odniesieniu do dopasowania wielomianów trygonometrycznych, i tak zazwyczaj radzę sobie z tego rodzaju danymi (chyba że odpowiedź nie jest bardzo płynna, w takim przypadku przełączam się na proces Gaussa z okresowym jądrem).

Splajny są używane w modelowaniu regresji do modelowania możliwie złożonych, nieliniowych form funkcjonalnych. Trend wygładzony splajnem składa się z częściowo ciągłych wielomianów, których wiodący współczynnik zmienia się w każdym punkcie przerwania lub węźle. Splajn można określić zarówno pod względem wielomianu trendu, jak i punktów przerwania. Reprezentacja splajnu współzmiennej rozciąga pojedynczy wektor obserwowanych wartości do macierzy, której wymiarem jest stopień wielomianowy plus liczba węzłów.

Okresowe wersja wypustów jest tylko okresowe wersja usprawiedliwienia: dane są cięte na powtórzeń o długości okresu. Na przykład modelowanie trendu dobowego w wielodniowym eksperymencie na szczurach wymagałoby przekodowywania czasu eksperymentu na przyrosty 24-godzinne, więc 154 godzina byłaby wartością modulo 24 równą 10 (154 = 6 * 24 + 10). Jeśli dopasujesz regresję liniową do danych cięcia, oszacuje ona przebieg fali piły dla trendu. Jeśli wpiszesz funkcję kroku gdzieś w tym okresie, będzie to kwadratowy przebieg pasujący do szeregu. Splajn jest w stanie wyrazić znacznie bardziej wyrafinowaną falkę. Za to, co jest warte, w splinespakiecie znajduje się funkcja, periodicSplinektóra właśnie to robi.

$p$$n_k$$p$$p+i$$i \le n_k$$S_{p+i} = (X - k_i)^p\mathcal{I}(X<k_i)$$k$

myspline <- function(x, degree, knots) {
  knots <- sort(knots)
  val <- cbind(x, outer(x, knots, `-`))
  val[val < 0] <- 0
  val <- val^degree
  if(degree > 1)
    val <- cbind(outer(x, 1:{degree-1}, `^`), val)
  colnames(val) <- c(
    paste0('spline', 1:{degree-1}, '.1'),
    paste0('spline', degree, '.', seq(length(knots)+1))
  )
  val
}

$2\pi$$\tau$

x <- seq(0, 2*pi, by=pi/2^8)
y <- sin(x)
plot(x,y, type='l')
s <- myspline(x, 2, pi)
fit <- lm(y ~ s)
yhat <- predict(fit)
lines(x,yhat)

Zobaczysz, że są dość zgodne. Ponadto konwencja nazewnictwa umożliwia interpretację. W wyniku regresji widać:

> summary(fit)

Call:
lm(formula = y ~ s)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.04564 -0.02050  0.00000  0.02050  0.04564 

Coefficients:
             Estimate Std. Error  t value Pr(>|t|)    
(Intercept) -0.033116   0.003978   -8.326 7.78e-16 ***
sspline1.1   1.268812   0.004456  284.721  < 2e-16 ***
sspline2.1  -0.400520   0.001031 -388.463  < 2e-16 ***
sspline2.2   0.801040   0.001931  414.878  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02422 on 509 degrees of freedom
Multiple R-squared:  0.9988,    Adjusted R-squared:  0.9988 
F-statistic: 1.453e+05 on 3 and 509 DF,  p-value: < 2.2e-16

$\pi/2$

Zakładam, że znasz częstotliwość dostępnych danych. Jeśli w danych brakuje komponentu wzrostu lub średniej ruchomej, możesz przekształcić długi szereg czasowy w repliki krótkich szeregów trwających 1 okres. Masz teraz repliki i możesz użyć analizy danych do oszacowania powtarzającego się trendu.

Załóżmy, że generuję następujące nieco głośne, bardzo długie szeregi czasowe:

x <- seq(1, 100, by=0.01)
y <- sin(x) + rnorm(length(x), 0, 10)
xp <- x %% (2*pi)
s <- myspline(xp, degree=2, knots=pi)
lm(y ~ s)

Wynikowa wydajność pokazuje rozsądną wydajność.

> summary(fit)

Call:
lm(formula = y ~ s)

Residuals:
    Min      1Q  Median      3Q     Max 
-39.585  -6.736   0.013   6.750  37.389 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.48266    0.38155  -1.265 0.205894    
sspline1.1   1.52798    0.42237   3.618 0.000299 ***
sspline2.1  -0.44380    0.09725  -4.564 5.09e-06 ***
sspline2.2   0.76553    0.18198   4.207 2.61e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.949 on 9897 degrees of freedom
Multiple R-squared:  0.006406,  Adjusted R-squared:  0.006105 
F-statistic: 21.27 on 3 and 9897 DF,  p-value: 9.959e-14