【问题标题】:interpretation of the output of R function bs() (B-spline basis matrix)R函数bs()输出的解释(B样条基矩阵)
【发布时间】:2019-06-19 06:04:18
【问题描述】:

我经常使用 B 样条进行回归。到目前为止,我从来不需要详细了解bs 的输出:我只需选择我感兴趣的模型,并将其与lm 匹配即可。但是,我现在需要在外部(非 R)代码中重现 b 样条模型。那么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

编辑:这与this question 相似但不完全相同。该问题询问回归系数的解释,但我比这更早一步:我想了解模型矩阵系数的含义。如果我尝试按照the first answer 中的建议制作相同的图,我会得到一个混乱的图:

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 样条基函数,因为它们的节点太多(每个函数应该只有一个)。

The second answer 实际上允许我在 R 之外重建我的模型,所以我想我可以这样做。但是,该答案也不能完全解释b 矩阵的元素是什么:它处理线性回归的系数,我还没有在这里介绍。这确实是我的最终目标,但我也想了解这个中间步骤。

【问题讨论】:

  • @ZheyuanLi,嗯,不。问题是关于lm 系数,我问的是基函数。 This answer 没有解释矩阵的单个系数是什么。如果我制作第一个答案中建议的相同类型的图,我会得到垃圾(绝对不是 B 样条函数,因为它们的最大值应该是 1)。 The other answer 更好,但仍不完全符合我的要求。我编辑了问题以说明原因。
  • 哦,等等!我知道了!我要回答我自己的问题,现在我明白那个矩阵是什么了:)
  • 首先,抱歉,我没有注意到您是答案的作者,否则我会这样称呼您。也许我们对基函数有不同的术语。对我来说,基函数是无限维对象(函数),你的答案显示了它们。但是矩阵b 的列不是基函数,而是样本点x &lt;- c(0.0, 11.0, 17.9, 49.3, 77.4) 中的基函数获得的值。 ctd..
  • ..ctd 我认为这并不完全相同:事实上,如果我在我的情况下重复你的情节,它就会变得一团糟(请参阅我编辑的问题)。但是,如果我不是针对向量 x 绘制 B 样条函数,而是针对节点向量绘制 B 样条函数,那么我会得到相同的图。在你的回答中,如果我解释正确,样本点和节点是相同的,所以这个问题不会出现。然而,在我的情况下,它们不是,这就是为什么我的矩阵b 包含不同于01 的元素。也许对你来说,这两个案例很明显是相同的,但我看不到。

标签: r matrix bspline


【解决方案1】:

矩阵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  

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

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

内部结是

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

那么三个基函数,如here所示,分别是:

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")

现在,考虑b[,1]

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

这些必须是x &lt;- 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] &gt; knots[3] = 38.83333b1 在那里为 0。其他两列可以给出类似的解释。

【讨论】:

  • 感谢您撰写本文!很有帮助
  • 您能解释一下为什么要用 x 值除以它们的结吗?你说“the value must be 11/13.3 = 0.8270677, as expected.”,但我不确定你为什么要划分这两个值。
  • @Parseltongue 这只是 x = 11 的值 (0, 0) 和 (13.3, 1) 的线性插值。在区间[0, 13.3]内,基函数b1是线性的,不是吗?
【解决方案2】:

只是对上面@DeltaIV 的出色答案的一个小修正(看起来我无法评论。)

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

注意b1 应该是这样的

\frac{t}{13.3}I(0&lt;=t&lt;13.3)+\frac{38.83333-t}{38.83333-13.3}I(13.3&lt;=t&lt;38.83333)

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2019-12-27
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2018-10-19
    • 2016-09-18
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多