週期樣條擬合週期數據
在對此問題的評論中，用戶@whuber 引用了使用周期版本的樣條曲線來擬合週期數據的可能性。我想更多地了解這種方法，特別是定義樣條的方程，以及如何在實踐中實現它們（我主要是
R用戶，但如果需要，我可以使用 MATLAB 或 Python）。此外，但這是一個“很高興擁有”，很高興了解三角多項式擬合的可能優點/缺點，這是我通常處理這類數據的方式（除非響應不是很平滑，在這種情況下，我切換到具有周期性內核的高斯過程）。
樣條曲線用於回歸建模以對可能複雜的非線性函數形式進行建模。樣條平滑趨勢由分段連續多項式組成，其前導係數在每個斷點或節點處發生變化。可以根據趨勢的多項式次數以及斷點來指定樣條。協變量的樣條表示將觀察值的單個向量擴展為一個矩陣，該矩陣的維數是多項式次數加上節點數。
樣條曲線的周期性版本僅僅是任何回歸的周期性版本：數據被切割成周期長度的重複。因此，例如，在大鼠的多日實驗中模擬晝夜趨勢需要將實驗時間重新編碼為 24 小時增量，因此第 154 小時將是 10 的模 24 值（154 = 6*24 + 10）。如果您對切割數據進行線性回歸，它將估計趨勢的鋸齒波形。如果您在周期的某處擬合階躍函數，它將是適合該系列的方波。樣條能夠表達更複雜的小波。對於它的價值，在
splines包中，有一個功能periodicSpline正是這樣做的。我沒有發現 R 的默認樣條“bs”實現對解釋有用。所以我在下面寫了自己的腳本。對於度數樣條和結，這個表示給出了第一個列標準多項式表示，-th 列 () 被簡單地評估為在哪裡是節點的實際向量。
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 }對於一個小案例研究，在 0 的域上插入正弦趨勢到（或者) 像這樣：
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)你會發現他們很默契。此外，命名約定可以解釋。在回歸輸出中，您會看到：
> 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我的 spline1.1 度的第一組協變量是第一個斷點後面的第一個域的多項式趨勢。線性項是原點處切線的斜率，X=0。這幾乎是 1，這將由正弦曲線的導數 (cos(0) = 1) 表示，但我們必須記住，這些是近似值，並且外推二次趨勢的誤差容易出錯。二次項表示負的凹形。spline2.2 項表示與第一個二次斜率的差異，導致 0.4 正前導係數表示向上的凸形。所以我們現在有了可用於樣條輸出的解釋，並且可以相應地判斷推斷和估計。
我將假設您知道手頭數據的周期性。如果數據缺少增長或移動平均組件，您可以將較長的時間序列轉換為持續時間為 1 個期間的較短序列的複製品。您現在有了重複，可以使用數據分析來估計重複趨勢。
假設我生成了以下有點嘈雜、很長的時間序列：
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)結果輸出顯示了合理的性能。
> 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