smooth.spline()：拟合模型与用户指定的自由度不匹配

Question

这是我的代码运行

fun <- function(x) {1 + 3*sin(4*pi*x-pi)}
set.seed(1)
num.samples <- 1000
x <- runif(num.samples)
y <- fun(x) + rnorm(num.samples) * 1.5
fit <- smooth.spline(x, y, all.knots=TRUE, df=3)

尽管df=3，当我检查拟合模型时，输出是

Call:
smooth.spline(x = x, y = y, df = 3, all.knots = TRUE)
Smoothing Parameter  spar= 1.499954  lambda= 0.002508571 (26 iterations)
Equivalent Degrees of Freedom (Df): 9.86422

有人可以帮忙吗？谢谢！

Answer 1

请注意，从 R-3.4.0 (2017-04-21) 开始，smooth.spline 可以接受通过新添加的参数 lambda 直接指定 λ。但是在估算的时候还是会转换成内部的spar。所以下面的答案不受影响。

平滑参数λ / spar位于平滑控制的中心

平滑度由平滑参数 λ 控制。smooth.spline() 使用内部平滑参数 spar 而不是 λ:

spar = s0 + 0.0601 * log(λ)

为了进行不受约束的最小化，这种对数变换是必要的，例如GCV/CV。用户可以指定 spar 来间接指定 λ。当spar线性增长时，λ将呈指数增长。因此很少需要使用大 spar 值。

自由度df，也是根据λ:

定义的

其中 X 是基于 B-spline 的模型矩阵，S 是惩罚矩阵。

您可以检查它们与您的数据集的关系：

spar <- seq(1, 2.5, by = 0.1)
a <- sapply(spar, function (spar_i) unlist(smooth.spline(x, y, all.knots=TRUE, spar = spar_i)[c("df","lambda")]))

让我们画出 df ~ spar、λ ~ spar 和 log(λ) ~ spar:

par(mfrow = c(1,3))
plot(spar, a[1, ], type = "b", main = "df ~ spar",
     xlab = "spar", ylab = "df")
plot(spar, a[2, ], type = "b", main = "lambda ~ spar",
     xlab = "spar", ylab = "lambda")
plot(spar, log(a[2,]), type = "b", main = "log(lambda) ~ spar",
     xlab = "spar", ylab = "log(lambda)")

注意λ与[=19=的激进增长]，log(λ)与spar的线性关系，以及df与df之间相对平滑的关系spar.

smooth.spline() 拟合迭代 spar

如果我们手动指定 spar 的值，就像我们在 sapply() 中所做的那样，则不会为选择 spar 进行拟合迭代；否则 smooth.spline() 需要遍历多个 spar 值。如果我们

指定cv = TRUE / FALSE，拟合迭代旨在最小化CV/GCV分数；
指定df = mydf，拟合迭代旨在最小化(df(spar) - mydf) ^ 2。

最小化 GCV 很容易理解。我们不关心GCV分数，只关心对应的spar。相反，在最小化(df(spar) - mydf)^2时，我们往往关心迭代结束时的df值而不是spar！但请记住，这是一个最小化问题，我们永远无法保证最终的 df 与我们的目标值 mydf.

匹配

为什么你放了df = 3，却得到了df = 9.864?

迭代结束，可能意味着达到最小值，或达到搜索边界，或达到最大迭代次数。

我们离最大迭代次数限制还很远（默认500）；但我们没有达到最低要求。好吧，我们可能会到达边界。

不关注df，想想spar。

smooth.spline(x, y, all.knots=TRUE, df=3)$spar   # 1.4999

根据 ?smooth.spline，默认情况下，smooth.spline() 在 [-1.5, 1.5] 之间搜索 spar。即，当您放置 df = 3 时，最小化在搜索边界处终止，而不是达到 df = 3.

再次查看我们的 df 和 spar 之间的关系图。从图中可以看出，我们需要一些接近 2 的 spar 值才能得到 df = 3.

让我们使用control.spar参数：

fit <- smooth.spline(x, y, all.knots=TRUE, df=3, control.spar = list(high = 2.5))
# Smoothing Parameter  spar= 1.859066  lambda= 0.9855336 (14 iterations)
# Equivalent Degrees of Freedom (Df): 3.000305

现在你看，你最终得到 df = 3。我们需要一个 spar = 1.86.

更好的建议：不要使用all.knots = TRUE

看，你有1000条数据。使用 all.knots = TRUE 您将使用 1000 个参数。希望以 df = 3 结束意味着 1000 个参数中的 997 个被抑制。想象一下 λ 因此 spar 你需要多大！

尝试改用惩罚回归样条。把200个参数压成3个肯定轻松很多：

fit <- smooth.spline(x, y, nknots = 200, df=3)  ## using 200 knots
# Smoothing Parameter  spar= 1.317883  lambda= 0.9853648 (16 iterations)
# Equivalent Degrees of Freedom (Df): 3.000386

现在，您最终 df = 3 没有 spar 控制权。

smooth.spline()：拟合模型与用户指定的自由度不匹配

smooth.spline(): fitted model does not match user-specified degree of freedom

regression

r

spline

smoothing