r - `smooth.spline` 严重欠拟合长(周期)时间序列

Question

我想在 R 中平滑非常长、嘈杂的数据。但我发现对于高度周期性的数据，开箱即用的 smooth.spline() 很快就会崩溃，平滑后的数据开始出现振铃现象。

考虑一个余弦时间序列(有或没有噪声)

t <- seq(0,100*2*pi,length.out=3000)
y <- cos(t)# + rnorm(length(t), 0,0.05)

y100_s <- smooth.spline(y)$y

plot( y~t, type="l" )
lines( y100_s~t, col="blue" )

我们可以检查向 smooth.spline() 添加更多值的效果，

# rms increases as points are added to smooth.spline
rms <- sapply( seq(250,3000,by=250), function(i)
  sqrt( mean( (y[1:i] - smooth.spline(y[1:i])$y)^2 )) )

plot(rms)

即使在较低的频率下，拟合也会响起(可选)。

t <- seq(0,50*2*pi,length.out=3000)
y <- cos(t)# + rnorm(length(t), 0,0.05)

y50_s <- smooth.spline(y)$y

require(pracma)

peaks <- list(findpeaks(y50_s),findpeaks(-y50_s))

plot( y~t, type="l" )
lines( y50_s~t, col="red" )

lines( peaks[[1]][,1]~t[peaks[[1]][,2]], type="l" )
lines( -peaks[[2]][,1]~t[peaks[[2]][,2]], type="l" )

经过一番探索，此行为似乎是 spar 参数的函数，但我无法将其设置为足够小的值来消除这种影响。这可能是样条拟合的一个明显结果，也是依赖开箱即用方法的错误，但我希望能提供一些见解。是否有我可以在 smooth.spline() 中指定的控件，或用于平滑的替代建议/策略？

Answer 1

不知道你是不是一直在拟合周期信号。如果是这种情况，使用 mgcv::gam 中的周期性样条会好得多。但是，让我们暂时忘记这个问题。

如果您的数据具有高频率振荡，则必须选择足够数量的结，即合适的结密度，否则只会导致过度平滑(即欠拟合)。

看看你的例子:

t <- seq(0, 100 * 2 * pi, length.out = 3000)
y <- cos(t) # + rnorm(length(t), 0, 0.05)
fit <- smooth.spline(t, y)

您有 n = 3000 个数据点。默认情况下，当 n > 49 时，smooth.spline 使用的节点数比数据少得多。准确地说，它是由服务例程 .nknots.smspl 选择的。但这没有最优性的理由。所以这是否合理由你来论证。让我们检查一下:

length(fit$fit$nk) - 2L  ## or `.nknots.smspl(3000)`
# [1] 194

fit$df
# [1] 194

它仅使用 194 节，模型最终具有 194 个自由度，没有惩罚效果。正如我之前所说，你最终会出现欠拟合:

plot(t, y, type = "l", col = "gray")
lines(fit, col = 2)

理想情况下，惩罚回归最终的自由度远小于结数。人们常常忘记惩罚是用来修复由原始非惩罚回归导致的过拟合问题。如果我们连惩罚效果都看不到，那么原始的非惩罚模型是欠拟合数据，所以增加结数直到我们达到过拟合状态。如果你懒得想这个，设置all.knots = TRUE。单变量平滑样条的计算成本非常低，成本为 O(n)。即使您将所有数据都用作节点，您也不会遇到效率问题。

fit <- smooth.spline(t, y, all.knots = TRUE)

length(fit$fit$nk) - 2L
# [1] 3000

fit$df
# [1] 3000

哦，我们还是没有看到惩罚的效果，为什么呢？因为我们没有嘈杂的数据。你没有给你的 y 添加噪音，所以通过使用所有的结，我们正在做插值。向 y 添加一些噪音，以真正理解我对惩罚的解释。

set.seed(0)
t <- seq(0, 100 * 2 * pi, length.out = 3000)
y <- cos(t) + rnorm(length(t), 0, 0.05)

fit <- smooth.spline(t, y, all.knots = TRUE)

length(fit$fit$nk)
# [1] 3002

fit$df
# [1] 705.0414

请注意 705 与 3000 相比小了多少。看过拟合样条吗？

plot(t, y, type = "l", col = "gray")
lines(fit, col = 2)

既没有欠拟合也没有过拟合；惩罚导致偏差和方差之间的最佳权衡。