Periodische Splines zur Anpassung an periodische Daten

In einem Kommentar zu dieser Frage zitierte Benutzer @whuber die Möglichkeit, eine periodische Version von Splines zu verwenden, um periodische Daten anzupassen. Ich würde gerne mehr über diese Methode erfahren, insbesondere über die Gleichungen, die die Splines definieren, und wie man sie in der Praxis implementiert (ich bin meistens ein RBenutzer, kann aber bei Bedarf mit MATLAB oder Python auskommen). Auch, aber dies ist ein "nice to have", wäre es großartig, über mögliche Vor- und Nachteile in Bezug auf die Anpassung trigonometrischer Polynome Bescheid zu wissen, wie ich normalerweise mit dieser Art von Daten umgehe (es sei denn, die Antwort ist nicht sehr glatt). In diesem Fall wechsle ich zum Gaußschen Prozess mit periodischem Kernel.

Splines werden in der Regressionsmodellierung verwendet, um möglicherweise komplexe, nichtlineare Funktionsformen zu modellieren. Ein Spline-Glättungstrend besteht aus stückweise kontinuierlichen Polynomen, deren Leitkoeffizient sich an jedem Haltepunkt oder Knoten ändert. Der Spline kann sowohl hinsichtlich des Polynomgrades des Trends als auch der Haltepunkte angegeben werden. Eine Spline-Darstellung einer Kovariate erweitert einen einzelnen Vektor beobachteter Werte in eine Matrix, deren Dimension der Polynomgrad plus die Anzahl der Knoten ist.

Eine periodische Version von Splines ist lediglich eine periodische Version einer Regression: Die Daten werden in Replikate der Länge der Periode geschnitten. So würde beispielsweise die Modellierung eines Tagesverlaufs in einem mehrtägigen Experiment an Ratten eine Rekodierungszeit des Experiments in Schritten von 24 Stunden erfordern, sodass die 154. Stunde der Modulo-24-Wert von 10 (154 = 6 * 24 + 10) wäre. Wenn Sie eine lineare Regression an die Schnittdaten anpassen, wird eine Sägezahnwellenform für den Trend geschätzt. Wenn Sie eine Schrittfunktion irgendwo in der Periode anpassen, ist dies eine Rechteckwellenform, die zur Serie passt. Der Spline kann ein viel komplexeres Wavelet ausdrücken. Für das, was es wert ist splines, gibt es im Paket eine Funktion, periodicSplinedie genau dies tut.

Ich finde Rs Standard-Spline-Implementierung "bs" nicht nützlich für die Interpretation. Also habe ich unten mein eigenes Skript geschrieben. Für einen Spline vom Grad mit Knoten gibt diese Darstellung den ersten Spalten die Standardpolynomdarstellung, die ten Spalten ( ) werden einfach als ) bewertet wobei der tatsächliche Knotenvektor ist.n k p p + i i ≤ n k S p + i = ( X - k i ) p I ( X < k i ) $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
}

Interpolieren Sie für eine kleine Fallstudie einen sinusförmigen Trend in der Domäne von 0 bis (oder ) wie folgt:$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)

Sie werden sehen, dass sie ziemlich übereinstimmen. Ferner ermöglicht die Namenskonvention die Interpretation. In der Regressionsausgabe sehen Sie:

> 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

Der erste Satz von Kovariaten für meinen Spline1.1-Grad ist der Polynomtrend für die erste Domäne hinter dem ersten Haltepunkt. Der lineare Term ist die Steigung der Tangente am Ursprung, X = 0. Dies ist fast 1, was durch die Ableitung der Sinuskurve (cos (0) = 1) angezeigt würde, aber wir müssen bedenken, dass dies Näherungen sind und der Fehler der Extrapolation des quadratischen Trends out anfällig ist zum Fehler. Der quadratische Term zeigt eine negative, konkave Form an. Der Ausdruck spline2.2 gibt einen Unterschied zur ersten quadratischen Steigung an, der zu einem positiven Leitkoeffizienten von 0,4 führt, der eine nach oben gerichtete konvexe Form anzeigt. Wir haben jetzt also eine Interpretation für die Spline-Ausgabe und können die Inferenz und Schätzungen entsprechend beurteilen.$\pi/2$

Ich gehe davon aus, dass Sie die Periodizität der vorliegenden Daten kennen. Wenn den Daten eine Wachstums- oder gleitende Durchschnittskomponente fehlt, können Sie eine lange Zeitreihe in Replikate einer kurzen Reihe mit einer Dauer von 1 Periode umwandeln. Sie haben jetzt Replikate und können die Datenanalyse verwenden, um den wiederkehrenden Trend abzuschätzen.

Angenommen, ich generiere die folgenden etwas lauten, sehr langen Zeitreihen:

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)

Die resultierende Ausgabe zeigt eine angemessene Leistung.

> 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