在SciPy中生成BSpline基，如R中的bs（）

Question

对于N个一维数据X，我想计算K个三次B样条曲线上的每个点。在R中，有一个简单的函数，具有直观的API，称为bs。实际上有一个python包patsy，它replicates this，但是我不能使用这个包——只有scipy之类的

在浏览了关于样条曲线相关函数的scipy.interpolate文档之后，我能找到的最接近的是BSpline，或BSpline.basis_元素，但如何只获得K个基函数对我来说是完全神秘的。我尝试了以下方法：

import numpy as np
import scipy.interpolate as intrp
import matplotlib.pyplot as plt
import patsy # for comparison

# in Patsy/R: nice and sensible
x = np.linspace(0., 1., 100)
y = patsy.bs(x, knots=np.linspace(0,1,4), degree=3)
plt.subplot(1,2,1)
plt.plot(x,y)
plt.title('B-spline basis')

# in scipy: ?????
y_py = np.zeros((x.shape[0], 6))
for i in range(6):
    y_py[:,i] = intrp.BSpline(np.linspace(0,1,10),(np.arange(6)==i).astype(float), 3, extrapolate=False)(x)


plt.subplot(1,2,2)
plt.plot(x,y_py)
plt.title('Something else')

它不起作用，让我意识到我实际上不知道这个函数在做什么。首先，它不会接受少于8个内部结，我不明白为什么。其次，它只认为样条曲线定义在（1/3，2/3）ish范围内，这可能意味着它出于某种原因忽略了前3个和最后3个节点值？我需要打绳结吗

任何帮助都将不胜感激

编辑：我已经解决了这个矛盾，实际上，BSpline似乎忽略了结的前3个和最后3个值。我仍然想知道为什么会有这种差异，这样我就不会因为花了那么多时间调试一个奇怪的界面而感到难过

对于后人来说，以下是生成基函数的代码

import numpy as np
import scipy.interpolate as intrp
import matplotlib.pyplot as plt
import patsy # for comparison

these_knots = np.linspace(0,1,5)

# in Patsy/R: nice and sensible
x = np.linspace(0., 1., 100)
y = patsy.bs(x, knots=these_knots, degree=3)
plt.subplot(1,2,1)
plt.plot(x,y)
plt.title('B-spline basis')

# in scipy: ?????
numpyknots = np.concatenate(([0,0,0],these_knots,[1,1,1])) # because??
y_py = np.zeros((x.shape[0], len(these_knots)+2))
for i in range(len(these_knots)+2):
    y_py[:,i] = intrp.BSpline(numpyknots, (np.arange(len(these_knots)+2)==i).astype(float), 3, extrapolate=False)(x)

plt.subplot(1,2,2)
plt.plot(x,y_py)
plt.title('In SciPy')

Answer 1

看起来您已经找到了答案，但是为了澄清为什么需要在边缘定义多个节点，您可以阅读scipy docs。它们是使用Cox de Boor递归公式定义的。该公式首先定义给定结点之间的相邻支撑域，其常数值为1（零阶）。这些函数被卷积以获得高阶基函数。因此，两个域构成一个一阶基函数，三个域构成一个二阶基函数，四个域（=5个节点）构成一个三阶基函数，在这5个节点的范围内得到支持。如果想要阶数为k=3的n个基函数，则需要有（n+k+1）个节点

最小速度为8节时，n>；=k+1，它给出2*（k+1）。基本间隔t[k]。。。scipy中的t[n]是唯一可以定义全阶基函数的范围。为确保该基本间隔到达外部结点，两端结通常具有（k+1）的多重性。可能scipy只在“其他东西”结果中显示了这个基本间隔

请注意，您还可以使用

y_py[:,i] = intrp.BSpline.basis_element(numpyknots[i:i+5], extrapolate=False)(x)

这也消除了x=1时的差异