R函数bs()(B样条基矩阵)输出的解释

15

我经常使用 B-Splines 进行回归。到目前为止,我从未需要详细了解 bs 的输出:我只会选择我感兴趣的模型,并用 lm 拟合它。然而,我现在需要在外部(非 R )代码中复现一个B-spline模型。那么,bs 生成的矩阵意味着什么?例如:

x <- c(0.0, 11.0, 17.9, 49.3, 77.4)
bs(x, df = 3, degree = 1) # generate degree 1 (linear) B-splines with 2 internal knots
#              1         2         3
# [1,] 0.0000000 0.0000000 0.0000000    
# [2,] 0.8270677 0.0000000 0.0000000    
# [3,] 0.8198433 0.1801567 0.0000000    
# [4,] 0.0000000 0.7286085 0.2713915    
# [5,] 0.0000000 0.0000000 1.0000000   
# attr(,"degree")
# [1] 1
# attr(,"knots")
# 33.33333% 66.66667% 
#  13.30000  38.83333 
# attr(,"Boundary.knots")
# [1]  0.0 77.4
# attr(,"intercept")
# [1] FALSE
# attr(,"class")
# [1] "bs"     "basis"  "matrix"

好的,所以我在输入中指定了degree为1。 knots 告诉我两个内部结节分别位于 x = 13.3000 和 x = 38.8333 处,看到结节位于固定的分位数上有些惊讶。我希望 R 能够找到最适合我的数据的分位数,但这会使模型不再线性,并且在不知道响应数据的情况下是不可能的。intercept = FALSE 表示基函数中未包括截距(这是件好事吗?我一直被教导不要在没有截距的情况下拟合线性模型...好吧,猜想 lm 仍然会添加一个截距)。

然而,矩阵怎么办?我不太明白如何解释它。有三列,我认为意味着基函数有三个。这是有道理的:如果我有两个内部结节 K1K2 ,那么我将在左边界结节 B1K1 之间拥有一个样条,另一个样条在 K1K2 之间,最后一个样条在 K2B2 之间,所以......三个基函数,好的。但是哪些是基函数呢?例如,这一列表示什么?

#              1
# [1,] 0.0000000
# [2,] 0.8270677
# [3,] 0.8198433
# [4,] 0.0000000
# [5,] 0.0000000

编辑:这个问题与这个问题类似但不完全相同。那个问题问的是回归系数的解释,但我想先了解模型矩阵系数的含义。如果我尝试按照第一个答案建议的制作相同的图表,我会得到一个混乱的图表:

b <- bs(x, df = 3, degree = 1)
b1 <- b[, 1]  ## basis 1
b2 <- b[, 2]  ## basis 2
b3 <- b[,3]
par(mfrow = c(1, 3))
plot(x, b1, type = "l", main = "basis 1: b1")
plot(x, b2, type = "l", main = "basis 2: b2")
plot(x, b3, type = "l", main = "basis 3: b3")

在此输入图片描述

这些不能是B样条基函数,因为它们有太多的节点(每个函数应该只有一个节点)。

第二个答案实际上可以让我在R之外重建我的模型,所以我想我可以选择那个答案。然而,也就是那个答案并没有准确地解释b矩阵中的元素:它涉及线性回归的系数,在这里我还没有介绍过。确实,那是我的最终目标,但我也想理解这个中间步骤。


@ZheyuanLi,不是的。问题是关于“lm”系数的,我问的是基函数。这个答案没有解释矩阵的单个系数是什么。如果我按照第一个答案建议的方式绘制相同类型的图形,我会得到垃圾(肯定不是B样条函数,因为它们的最大值应该是1)。另一个答案更好,但仍然不完全是我所要求的。我编辑了问题以说明原因。 - DeltaIV
哦,等等!我明白了!我要回答自己的问题,现在我明白那个矩阵是什么意思了 :) - DeltaIV
首先,抱歉,我没有注意到您是答案的作者,否则我会称呼您的。也许我们对基函数有不同的术语。对我来说,基函数是无限维对象(函数),而您的答案展示了它们。但是,矩阵b的列在我看来不是基函数,而是基函数在样本点x <- c(0.0, 11.0, 17.9, 49.3, 77.4)中获得的值。继续... - DeltaIV
我认为这不完全相同:事实上,如果我在我的情况下重复您的图,它会变得混乱不堪(请参见我的编辑问题)。但是,如果我将B样条函数绘制到结节向量而不是向量x中,则可以获得您相同的绘图。在您的答案中,如果我理解正确,样本点和结节是相同的,因此这个问题不会出现。然而,在我的情况下,它们并不相同,这就是为什么我的矩阵b包含与01不同的元素的原因。也许对您来说,这两种情况是相同的,但我看不出来。 - DeltaIV
2个回答

12

矩阵b

#              1         2         3
# [1,] 0.0000000 0.0000000 0.0000000    
# [2,] 0.8270677 0.0000000 0.0000000    
# [3,] 0.8198433 0.1801567 0.0000000    
# [4,] 0.0000000 0.7286085 0.2713915    
# [5,] 0.0000000 0.0000000 1.0000000  

实际上,它只是每个点的三个基函数值的矩阵,这一点对我来说应该很明显,因为它与多项式线性模型的解释完全相同。事实上,由于边界结点是

bknots <- attr(b,"Boundary.knots")
# [1]  0.0 77.4

内部节点是

iknots <- attr(b,"knots")
# 33.33333% 66.66667% 
#  13.30000  38.83333 

然后,如这里所示,三个基本函数是:

knots <- c(bknots[1],iknots,bknots[2])
y1 <- c(0,1,0,0)
y2 <- c(0,0,1,0)
y3 <- c(0,0,0,1)
par(mfrow = c(1, 3))
plot(knots, y1, type = "l", main = "basis 1: b1")
plot(knots, y2, type = "l", main = "basis 2: b2")
plot(knots, b3, type = "l", main = "basis 3: b3")

enter image description here

现在,请考虑 b[,1]

#              1
# [1,] 0.0000000
# [2,] 0.8270677
# [3,] 0.8198433
# [4,] 0.0000000
# [5,] 0.0000000
这些必须是x <- c(0.0, 11.0, 17.9, 49.3, 77.4)b1的值。实际上,b1knots[1] = 0处为0,在knots[2] = 13.3000处为1,这意味着在x[2](11.0)中该值必须是11/13.3 = 0.8270677,正如预期的那样。同样,由于b1对于knots[3] = 38.83333为0,因此x[3](17.9)中的值必须是(38.83333-13.3)/17.9 = 0.8198433。由于x[4], x[5] > knots[3] = 38.83333,所以那里的b1为0。其他两列可以给出类似的解释。

3
谢谢您写下这篇文章!非常有帮助。 - Parseltongue
你能解释一下为什么要将x值除以它们的节点吗?你说:“该值必须为11/13.3 = 0.8270677,正如预期的那样。”,但我不确定为什么要将这两个值相除。 - Parseltongue
3
@Parseltongue,这只是在x=11时,对值(0,0)和(13.3,1)进行线性插值。在[0,13.3]区间内,基函数“b1”是线性的,是吗? - DeltaIV

3

对@DeltaIV上面的优秀回答进行小小的更正(看起来我不能发表评论)。

所以在b1中,当他计算b1(x [3])时,应该通过线性插值为(38.83333-17.9)/(38.83333-13.3)=0.8198433。其他都是完美的。

注意,在这种情况下,b1应该像这样:

\frac{t}{13.3}I(0<=t<13.3)+\frac{38.83333-t}{38.83333-13.3}I(13.3<=t<38.83333)


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接