从 R 中的 smooth.spline 检索/重现设计矩阵

Question

我想知道如何检索或复制在幕后使用的设计矩阵smooth.spline。似乎没有smooth.spline办法返回这个矩阵。因此，我尝试使用该bs函数来重现矩阵。

当我运行smooth.spline并通过参数指定结时all.knots，与通过将相同结传递给函数而获得的 bspline 矩阵相比，拟合的参数数量减少了 1 bs。smooth.spline有人可以帮助我理解如何通过给定一组已知的结来重现在引擎盖下使用的训练矩阵吗？谢谢！

下面是重现该行为的代码：

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

Answer 1

这并没有解释如何使用 的输出的内部数据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 ) = 西.

请注意，这样的 X 不是唯一的，因为如果我们将与 b 正交的向量添加到 X 的行中，则方程仍然满足。

虽然我们显示的 X 确实满足所要求的方程，但我不确定它是否真的是您想要的。请注意，查看源代码，我们可以将 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))
}
lM <- aux2Mat(mod$auxM)
A <- lM$XWX + mod$lambda * lM$Sigma
R <- Matrix::chol(A)

以下两行给出相同的结果：

bb <- solve(A, lM$XWy)
b <- mod$fit$coef

max(abs(b - bb))
## [1] 1.125863e-10

这两个也是一样的

crossprod(b, as.matrix(lM$XWX)) %*% b
##            [,1]
## [1,] 0.01597881
crossprod(X %*% b)
##            [,1]
## [1,] 0.01597881

这两个是一样的

y %*% X %*% b
##            [,1]
## [1,] 0.02039534
lM$XWy %*% b
##            [,1]
## [1,] 0.02039534

更新

更正。