r语言 - 在Python中获得' bs '(样条)等价

Question

在r编程语言中，下面的

require(stats)
require(splines)
knots = quantile(women$height, seq(0.1,0.9,length.out = 5))
bs(women$height, knots=knots, degree=3)

返回

1   2   3   4   5   6   7   8
0.0000000000    0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.0000000000    0.00000000
0.6284418529    0.323939099 0.024295432 0.000000000 0.000000000 0.000000000 0.0000000000    0.00000000
0.2155814707    0.599894720 0.182883868 0.001639942 0.000000000 0.000000000 0.0000000000    0.00000000
0.0349854227    0.495626822 0.438289602 0.031098154 0.000000000 0.000000000 0.0000000000    0.00000000
0.0001619695    0.245586330 0.620809038 0.133442663 0.000000000 0.000000000 0.0000000000    0.00000000
0.0000000000    0.072886297 0.584548105 0.338678328 0.003887269 0.000000000 0.0000000000    0.00000000
0.0000000000    0.009110787 0.384718173 0.561892614 0.044278426 0.000000000 0.0000000000    0.00000000
0.0000000000    0.000000000 0.166666667 0.666666667 0.166666667 0.000000000 0.0000000000    0.00000000
0.0000000000    0.000000000 0.044278426 0.561892614 0.384718173 0.009110787 0.0000000000    0.00000000
0.0000000000    0.000000000 0.003887269 0.338678328 0.584548105 0.072886297 0.0000000000    0.00000000
0.0000000000    0.000000000 0.000000000 0.133442663 0.620809038 0.245586330 0.0001619695    0.00000000
0.0000000000    0.000000000 0.000000000 0.031098154 0.438289602 0.495626822 0.0349854227    0.00000000
0.0000000000    0.000000000 0.000000000 0.001639942 0.182883868 0.599894720 0.2155814707    0.00000000
0.0000000000    0.000000000 0.000000000 0.000000000 0.024295432 0.323939099 0.6284418529    0.02332362
0.0000000000    0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.0000000000    1.00000000

是否有对应的Python ?我已经从scipy中尝试了BSpline，但它要求系数已经已知并传入。

如何生成b样条基矩阵，传入数组、结点和度数?

---

要复制输入的Python，可以这样做:

import numpy as np
women_height = np.array([58,59,60,61,62,63,64,65,66,67,68,69,70,71,72])
knots = array([59.4, 62.2, 65. , 67.8, 70.6])

Answer 1

将评论转换为答案，BSpline.design_matrix正在构建您所追求的内容，以csr稀疏格式。它将在scipy 1.8发布时提供。在此之前，您可以获取scipy的主分支，或者使用文档建议的解决方案(https://scipy.github.io/devdocs/reference/generated/scipy.interpolate.BSpline.design_matrix.html#scipy.interpolate.BSpline.design_matrix):

t = ... 
c = np.eye(len(t) - k - 1)
design_matrix_gh = BSpline(t, c, k)(x)

EDIT: R文档，https://www.rdocumentation.org/packages/splines/versions/3.6.2/topics/bs，声明knots参数是内部节。scipy的BSpline不会自动填充结，所以你需要自己做。使用OP数据:

In [22]: women_height = np.array([58,59,60,61,62,63,64,65,66,67,68,69,70,71,72])
...: knots = np.array([59.4, 62.2, 65. , 67.8, 70.6])
In [23]: t = np.r_[(58,)*4, knots, (72,)*4]   # <<<<<< here
In [24]: m = BSpline.design_matrix(women_height, t, k=3)
In [25]: with np.printoptions(linewidth=120, precision=5):
...:     print(m.toarray())
...:
[[1.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00]
[2.33236e-02 6.28442e-01 3.23939e-01 2.42954e-02 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00]
[0.00000e+00 2.15581e-01 5.99895e-01 1.82884e-01 1.63994e-03 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00]
[0.00000e+00 3.49854e-02 4.95627e-01 4.38290e-01 3.10982e-02 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00]
[0.00000e+00 1.61970e-04 2.45586e-01 6.20809e-01 1.33443e-01 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00]
[0.00000e+00 0.00000e+00 7.28863e-02 5.84548e-01 3.38678e-01 3.88727e-03 0.00000e+00 0.00000e+00 0.00000e+00]
[0.00000e+00 0.00000e+00 9.11079e-03 3.84718e-01 5.61893e-01 4.42784e-02 0.00000e+00 0.00000e+00 0.00000e+00]
[0.00000e+00 0.00000e+00 0.00000e+00 1.66667e-01 6.66667e-01 1.66667e-01 0.00000e+00 0.00000e+00 0.00000e+00]
[0.00000e+00 0.00000e+00 0.00000e+00 4.42784e-02 5.61893e-01 3.84718e-01 9.11079e-03 0.00000e+00 0.00000e+00]
[0.00000e+00 0.00000e+00 0.00000e+00 3.88727e-03 3.38678e-01 5.84548e-01 7.28863e-02 0.00000e+00 0.00000e+00]
[0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 1.33443e-01 6.20809e-01 2.45586e-01 1.61970e-04 0.00000e+00]
[0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 3.10982e-02 4.38290e-01 4.95627e-01 3.49854e-02 0.00000e+00]
[0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 1.63994e-03 1.82884e-01 5.99895e-01 2.15581e-01 0.00000e+00]
[0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 2.42954e-02 3.23939e-01 6.28442e-01 2.33236e-02]
[0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 1.00000e+00]]

看起来类似于R从OP对第一列取模得到的输出。从R的文档中，目前还不清楚它是如何填充结向量的，但如果你想要相同的输出，你可以把第一列切掉(m.toarray()[1:, :]或一些这样的)。