r - 平滑样条(): fitted model does not match user-specified degree of freedom

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值(value)。

自由度 df , 也定义为 λ :



哪里 X是具有 B 样条基和 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)")

请注意 λ 的急剧增长与 spar ， log(λ)之间的线性关系和 spar ，以及 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()搜索 spar之间 [-1.5, 1.5] .即，当您输入 df = 3 时，最小化终止于搜索边界，而不是命中 df = 3 .

看看我们的 df 之间的关系图和 spar ， 再次。从图中看，我们需要一些 spar值接近 2 以导致 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控制。