我想知道如何检索或复制smooth.spline年前引擎盖下使用的设计矩阵.smooth.spline似乎没有办法返回这个矩阵.因此，我正在try 使用bs函数来重新生成矩阵.
当我运行smooth.spline并通过all.knots参数指定 node 时，与通过将相同的 node 传递给bs函数而获得的B样条阵相比，来自Fit的参数的数量减少了1.有没有人能帮我理解如何重现smooth.spline个已知 node 所使用的训练矩阵？谢谢!
以下是重现该行为的代码:
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
这并没有解释如何使用smooth.spline的输出的内部数据来获得更平滑的矩阵，但我们可以通过将单位矩阵通过smooth.spline来获得该矩阵(这也是求线性运算符的矩阵的一般方法).另请参阅https://www.stat.cmu.edu/~cshalizi/dst/18/hw/02/smoother.matrix.R
set.seed(123)
n <- 200
y <- rnorm(n)
x <- runif(n)
mod <- smooth.spline(x, y, keep.stuff = TRUE)
S <- apply(diag(n), 2, function(y) fitted(smooth.spline(x, y, df = mod$df)))
max(abs(S %*% y - fitted(mod)))
## [1] 2.183287e-08
如果你只需要更平滑的矩阵的对角线，那么hatvalues(mod)就可以了.
max(abs(diag(S) - hatvalues(mod)))
## [1] 1.676781e-08
这将计算一个满足xB=Sy的矩阵X.
b <- mod$fit$coef
X <- tcrossprod(fitted(mod), b/c(crossprod(b)))
max(abs(X %*% b - S %*% y))
## [1] 4.959864e-09
方程成立是因为如果X=(Sy)b‘/(b’b)，则xb=((Sy)b‘/(b’b))b=(Sy)(b‘b)/(b’b)=Sy.
请注意，这样的X不是唯一的，因为如果我们将与b垂直的向量加到X的行上，那么等式仍然成立.
尽管我们展示的X确实满足所要求的等式，但我不确定它是否真的是您想要的.请注意，查看source，我们可以将mod$AUXM转换为如下所示的矩数组表:
aux2Mat <- function(auxM) {
stopifnot(is.list(auxM),
identical(vapply(auxM, class, ""),
setNames(rep("numeric", 4), c("XWy", "XWX", "Sigma", "R"))))
## requireNamespace("Matrix")# want sparse matrices
nk <- length(XWy <- auxM[["XWy"]])
list(XWy = XWy,
XWX = Matrix::bandSparse(nk, k= 0:3, diagonals= matrix(auxM[[ "XWX" ]], nk,4), symmetric=TRUE),
Sigma= Matrix::bandSparse(nk, k= 0:3, diagonals= matrix(auxM[["Sigma"]], nk,4), symmetric=TRUE))
}
L <- aux2Mat(mod$auxM)
此外，除了可能的尺寸外，所有4个参数都会产生相同的结果:
bb <- solve(L$XWX + mod$lambda * L$Sigma, L$XWy); bb
(L$XWX + mod$lambda * L$Sigma) %*% bb
L$XWy
mod$fit$coef