Splines periódicos para ajustar dados periódicos

Em um comentário a esta pergunta , o usuário @whuber citou a possibilidade de usar uma versão periódica de splines para ajustar dados periódicos. Gostaria de saber mais sobre esse método, em particular as equações que definem os splines e como implementá-los na prática (sou principalmente um Rusuário, mas posso me contentar com MATLAB ou Python, se necessário). Além disso, mas é um "bom ter", seria ótimo saber sobre possíveis vantagens / desvantagens em relação ao ajuste de polinômios trigonométricos, que é como eu geralmente lido com esse tipo de dados (a menos que a resposta não seja muito suave, nesse caso, mudo para o Processo Gaussiano com kernel periódico).

As splines são usadas na modelagem de regressão para modelar formas funcionais não lineares possivelmente complexas. Uma tendência suavizada de spline consiste em polinômios contínuos por partes cujo coeficiente principal muda em cada ponto de interrupção ou nó. O spline pode ser especificado em termos do grau polinomial da tendência, bem como dos pontos de interrupção. Uma representação spline de uma covariável estende um único vetor de valores observados para uma matriz cuja dimensão é o grau polinomial mais o número de nós.

Uma versão periódica de splines é apenas uma versão periódica de qualquer regressão: os dados são cortados em réplicas da duração do período. Assim, por exemplo, modelar uma tendência diurna em um experimento de sexta-feira em ratos exigiria o tempo de recodificação do experimento em incrementos de 24 horas; portanto, a 154ª hora seria o valor do módulo 24 de 10 (154 = 6 * 24 + 10). Se você ajustar uma regressão linear nos dados de corte, estimaria uma forma de onda dente de serra para a tendência. Se você ajustar uma função step em algum lugar do período, seria uma forma de onda quadrada que se encaixa na série. O spline é capaz de expressar uma wavelet muito mais sofisticada. Para o que vale a pena, no splinespacote, há uma função periodicSplineque faz exatamente isso.

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

Você verá que eles são bastante concordantes. Além disso, a convenção de nomenclatura permite a interpretação. Na saída de regressão, você vê:

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

Eu vou assumir que você conhece a periodicidade dos dados em mãos. Se os dados não apresentarem um componente de crescimento ou média móvel, você poderá transformar uma longa série temporal em réplicas de uma curta série com duração de 1 período. Agora você tem réplicas e pode usar a análise de dados para estimar a tendência recorrente.

Suponha que eu gere as seguintes séries temporais um tanto ruidosas e muito longas:

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)

A saída resultante mostra um desempenho razoável.

> 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