scipy - scipy.interpolate.InterpolatedUnivariateSpline.get_coeffs返回什么？

我尝试了以下方法：

spline= interpolate.InterpolatedUnivariateSpline(X, Y, k=3)
coefs= spline.get_coeffs()

在X和Y分别具有五个值的情况下，我最终得到了coefs
五个值。假设有五个数据点意味着四个样条曲线部分，并且
三次多项式有四个系数，我本来希望得到
四乘四= 16个系数。有谁知道如何解释get_coeffs方法返回的值？是否有记录在案的地方？

最佳答案

这些不是x，x ** 2等的系数：这样的单项式不适合表示样条曲线。相反，它们是B-splines的系数，这些系数是针对进行插值的特定网格计算的。所涉及的B样条曲线的数量等于数据点的数量，系数的数量也等于。例如，假设我们要对这些数据进行插值：

xv = [0, 1, 2, 3, 4, 5]
yv = [3, 6, 5, 7, 9, 1]

从度数k = 1（分段线性样条曲线）的简单情况开始。那么B样条就是这些“三角帽”功能：



有6个。它们中的每个在“其”网格点等于1，在所有其他网格点等于0。这使得编写内插样条线非常容易：它是y[0]*b[0] + ... + y[5]*b[5]。实际上，get_coeffs显示系数是y值本身。

InterpolatedUnivariateSpline(xv, yv, k=1).get_coeffs()  # [ 3.,  6.,  5.,  7.,  9.,  1.]

三次样条

现在它变得毛茸茸了，因为我们需要光滑的“帽子”，而不是上面的帽子。平滑度要求迫使它们更宽，因此每个B样条在几个网格点上都具有非零值。 （技术性：三次B样条曲线在最多3节的情况下具有非零值，但是在下图上，1和4尽管是网格点，但由于所谓的“非结”情况而并非结。 。）这是网格的B样条曲线：



为了获得这些，我使用了scipy.interpolate的旧方法splrep和splev，它们在后台调用了相同的fitpack例程。这里的系数向量是元组tck的第二项；我将其修改为一个1，其余0，从而创建基础样条线（b-spline）。

k = 3
tck = splrep(xv, yv, s=0, k=k)
xx = np.linspace(min(xv), max(xv), 500)
bsplines = []
for j in range(len(xv)):
    tck_mod = (tck[0], np.arange(len(xv)+2*k-2) == j, k)
    bsplines.append(splev(xx, tck_mod))
    plt.plot(xx, bsplines[-1])

现在我们有了一个列表bsplines，我们可以使用get_coeffs返回的系数将它们自己放在一起作为插值样条线：

coeffs = InterpolatedUnivariateSpline(xv, yv, k=3).get_coeffs()
interp_spline = sum([coeff*bspline for coeff, bspline in zip(coeffs, bsplines)])
plt.plot(xx, interp_spline)

如果您需要这些B样条曲线的公式，则B-splines上的Cox-de Boor递归公式可以提供帮助，但这是手工计算的繁琐工作。

SymPy可以给出formulas for B-splines，但是有一点点变化。一个人应该通过重复结束值来传递一组填充的结

[0, 0, 0, 0, 2, 3, 5, 5, 5, 5]

这是因为在0和5时，所有四个系数都改变值，而在1和4时，它们都没有改变，因此将其省略（“不打结”）。 （此外，SymPy（1.1.1）的当前版本存在重复打结的问题，但此问题将在下一版本中修复。）

from sympy import symbols, bspline_basis_set, plot
x = symbols('x')
xv_padded = [0, 0, 0, 0, 2, 3, 5, 5, 5, 5]
bs = bspline_basis_set(3, xv_padded, x)

现在，bs是一系列看上去很吓人的分段公式：

[Piecewise((-x**3/8 + 3*x**2/4 - 3*x/2 + 1, (x >= 0) & (x <= 2)), (0, True)),
Piecewise((19*x**3/72 - 5*x**2/4 + 3*x/2, (x >= 0) & (x <= 2)), (-x**3/9 + x**2 - 3*x + 3, (x >= 2) & (x <= 3)), (0, True)),
Piecewise((-31*x**3/180 + x**2/2, (x >= 0) & (x <= 2)), (11*x**3/45 - 2*x**2 + 5*x - 10/3, (x >= 2) & (x <= 3)), (-x**3/30 + x**2/2 - 5*x/2 + 25/6, (x >= 3) & (x <= 5)), (0, True)),
Piecewise((x**3/30, (x >= 0) & (x <= 2)), (-11*x**3/45 + 5*x**2/3 - 10*x/3 + 20/9, (x >= 2) & (x <= 3)), (31*x**3/180 - 25*x**2/12 + 95*x/12 - 325/36, (x >= 3) & (x <= 5)), (0, True)),
Piecewise((x**3/9 - 2*x**2/3 + 4*x/3 - 8/9, (x >= 2) & (x <= 3)), (-19*x**3/72 + 65*x**2/24 - 211*x/24 + 665/72, (x >= 3) & (x <= 5)), (0, True)),
Piecewise((x**3/8 - 9*x**2/8 + 27*x/8 - 27/8, (x >= 3) & (x <= 5)), (0, True))]