来自 smooth.spline 的更平滑矩阵

机器算法验证 r 平滑 样条
2022-04-12 12:18:09

我正在尝试获得不smooth.spline适合我的更平滑的矩阵。我使用来自http://statweb.stanford.edu/~tibs/ElemStatLearn/的骨矿物质密度数据。

bone <-read.table("bone.data", header=TRUE)
bmd_age <- smooth.spline(bone$age, bone$spnbmd, all.knots=TRUE, cv=TRUE)
bmd_fit <- predict(bmd_age, sort(bone$age))
df <- bmd_age$df

为了获得更平滑矩阵的列，我可以将响应向量 (bone$spnbmd) 替换为具有单个 1 的向量，其余的填充为 0。这既是教授推荐的，也是我在网上找到的https://stat.ethz.ch/pipermail/r-help/2006-June/108471.html。

所以我用

smooth.matrix = function(x){
  n = length(x);
  sm = matrix(0, n, n);  
  for(i in 1:n){
    y = rep(0, n); y[i]=1;
    sm_i = predict(smooth.spline(x, y, df=df),x)$y;
    sm[,i]= sm_i;
  }
  return(sm)
}

sm <- smooth.matrix(bone$age)

如果更平滑的矩阵是正确的，则以下两个量应该相同（均来自平滑样条模型的拟合值）。

fromsm <- sm%*%(bone$spnbmd[order(bone$age)])
fromfit <- bmd_fit$y

但是，它们不是。我认为问题在于smooth.matrix函数的定义，其中

sm_i = predict(smooth.spline(x, y, df=df),x)$y;

没有使用与 bmd_age 中相同的平滑拟合。我尝试过修复自由度、晶石、lambda、cv=FALSE 等，但到目前为止还没有运气。如何解决？

3个回答

经过几个小时的探索，我发现了以下内容：

因为 smooth.spline 算法选择 spar 而不是 lambda，所以只能（某种程度上）修复 spar。但是，lambda 是 spar 和另一个变量矩阵的函数。所以修复晶石不一定能修复 lambda。我还没有找到一种从smooth.spline 中提取更平滑矩阵的简单方法。但是，出于计算方差的目的，https ://stat.ethz.ch/pipermail/r-help/2006-June/108471.html（修复 spar 而不是 df）中提供的算法是对真实情况的近似估计更平滑的矩阵。从计算的方差，其中是估计的平滑矩阵，非常接近从正确的平滑矩阵计算的方差。SySTS

另一个 R 包“assist”有一个函数“ssr()”，它也可以进行平滑样条回归。它不如 smooth.spline 强大。但是内置函数“hat.ssr()”给出了从“ssr()”获得的模型的真正平滑矩阵。

上述答案仅近似平滑矩阵。这是一个解决方案，可以从 r 函数 smooth.spline() 获得精确的平滑矩阵。关键是要认识到平滑矩阵只是x和惩罚参数λ，允许我们平滑一个向量y~=(0,0,...,0,1,0,...,0)T，因此得到平滑矩阵的每一列。

library(splines)
x = seq(0, 100, by=0.1)
y = x*sin(x) + rnorm(length(x), 0, 0.1)

#use cross-validation to choose best smoothing parameter
spar = seq(0.01, 1, by = 0.01)
cv = rep_len(NA, length(spar))
for(i in 1:length(spar)){
    tempfit = smooth.spline(x, y, spar = spar[i], cv=TRUE, all.knots = TRUE)
    cv[i] = tempfit$cv.crit
}

#use the optimal smoothing parameter to produce a final fit
fit = smooth.spline(x, y, spar = spar[which(cv == min(cv))], cv=TRUE, all.knots = TRUE)

#calculate the smoothing matrix
L = matrix(nrow = length(x), ncol = length(x))
for(j in 1:length(x)){
    yi = rep_len(0, length(x))
    yi[j] = 1
    L[,j] = predict(smooth.spline(x, yi, lambda = fit$lambda, cv=TRUE,
                                  all.knots = TRUE), x)$y
}

矩阵L是得到的平滑矩阵。

接受的答案在这里不正确 -smooth.matrix工作得很好。

唯一的原因fromsm和fromfit上面的例子不一样是因为括号放错了。

替换fromsm <- sm%*%(bone$spnbmd[order(bone$age)])为fromsm <- (sm%*%bone$spnbmd)[order(bone$age)]，它们是相同的。

其它你可能感兴趣的问题
上一篇如何定义队列研究？ 下一篇1 对 1 线图的正式名称是什么？