Splines periódicos para ajustar datos periódicos

En un comentario a esta pregunta , el usuario @whuber citó la posibilidad de usar una versión periódica de splines para ajustar los datos periódicos. Me gustaría saber más acerca de este método, en particular las ecuaciones que definen las splines y cómo implementarlas en la práctica (soy principalmente un Rusuario, pero puedo hacerlo con MATLAB o Python, si surge la necesidad). Además, pero es "bueno tenerlo", sería bueno saber acerca de las posibles ventajas / desventajas con respecto al ajuste de polinomios trigonométricos, que es la forma en que generalmente trato con este tipo de datos (a menos que la respuesta no sea muy fluida, en cuyo caso cambio al Proceso Gaussiano con kernel periódico).

Las splines se usan en modelos de regresión para modelar formas funcionales posiblemente no complejas y complejas. Una tendencia suavizada de spline consiste en polinomios continuos por partes cuyo coeficiente principal cambia en cada punto de ruptura o nudo. La spline se puede especificar en términos del grado polinómico de la tendencia, así como los puntos de corte. Una representación spline de una covariable extiende un solo vector de valores observados en una matriz cuya dimensión es el grado polinómico más el número de nudos.

Una versión periódica de splines es simplemente una versión periódica de cualquier regresión: los datos se cortan en réplicas de la duración del período. Entonces, por ejemplo, modelar una tendencia diurna en un experimento de varios días en ratas requeriría un tiempo de grabación del experimento en incrementos de 24 horas, por lo que la hora 154 sería el valor del módulo 24 de 10 (154 = 6 * 24 + 10). Si ajusta una regresión lineal en los datos de corte, estimaría una forma de onda de diente de sierra para la tendencia. Si ajusta una función de paso en algún punto del período, sería una forma de onda cuadrada que se ajuste a la serie. La spline es capaz de expresar una wavelet mucho más sofisticada. Por lo que vale, en el splinespaquete, hay una función periodicSplineque hace exactamente esto.

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

Verás que son bastante concordantes. Además, la convención de nomenclatura permite la interpretación. En la salida de regresión ves:

> 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$

Asumiré que conoce la periodicidad de los datos disponibles. Si los datos carecen de un componente de crecimiento o promedio móvil, puede transformar una serie de tiempo larga en réplicas de una serie corta de una duración de 1 período. Ahora tiene réplicas y puede usar el análisis de datos para estimar la tendencia recurrente.

Supongamos que genero las siguientes series de tiempo algo ruidosas y muy largas:

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)

El resultado resultante muestra un rendimiento razonable.

> 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