如何获取/复现R中smooth.spline底层使用的设计矩阵

想了解如何获取或复现R中smooth.spline函数底层使用的设计矩阵，该函数没有返回此矩阵的途径。尝试用bs函数复现，但指定相同knots时，smooth.spline拟合的参数数量比bs生成的B样条矩阵列数少1，具体复现代码如下：

n = 200
y = rnorm(n=n)
x = runif(n=n)
knots = c(0,seq(0.001,0.999,.01),1)
mod = smooth.spline(x,y,all.knots=knots,keep.stuff=TRUE)

X = bs(x,knots=knots)

dim(X)[2] #105 
length(mod$fit$coef) #104

为什么参数数量差1？ #

smooth.spline默认使用的是自然三次样条，而splines包中的bs函数默认生成的是普通三次B样条。自然三次样条会在两个边界点（你的例子里是0和1）施加二阶导数为0的约束，这会减少1个自由度，对应参数数量就少1个。

如何复现smooth.spline的设计矩阵？ #

直接用splines包中的ns函数（自然样条函数）即可，它生成的就是自然三次样条的设计矩阵，列数会和smooth.spline的参数数量完全匹配：

library(splines)
# ns函数不需要把边界点（0和1）传入knots参数，手动指定Boundary.knots更清晰
X_ns = ns(x, knots = seq(0.001,0.999,.01), Boundary.knots = c(0,1))
dim(X_ns)[2] # 结果为104，和length(mod$fit$coef)一致

如果一定要用bs函数手动构造，需要额外施加自然样条的约束，但这种方法繁琐且容易出错，更推荐直接用ns函数。

另外补充：当你给smooth.spline指定all.knots时，它本质是做带惩罚的自然样条拟合；如果设置spar=0让惩罚项趋近于0，拟合出的系数就是无惩罚自然样条的系数，此时ns生成的矩阵就是对应的训练矩阵。

内容的提问来源于stack exchange，提问作者curious_dan