从节点和系数创建BSpline

Question

如果只知道点和系数，怎么能创建样条呢？我在这里使用scipy.interpolate.BSpline，但也对其他标准包开放。所以基本上，我想给一个人一个简单的系数数组，让他们能够重新建立数据的拟合。请看下面失败的红色虚线。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import BSpline, LSQUnivariateSpline

x = np.linspace(0, 10, 50) # x-data
y = np.exp(-(x-5)**2/4)    # y-data

# define the knot positions

t = [1, 2, 4, 5, 6, 8, 9]

# get spline fit

s1 = LSQUnivariateSpline(x, y, t)

x2 = np.linspace(0, 10, 200) # new x-grid
y2 = s1(x2) # evaluate spline on that new grid

# FAILED: try to construct BSpline using the knots and coefficients

k = s1.get_knots()
c = s1.get_coeffs()
s2 = BSpline(t,c,2)

# plotting

plt.plot(x, y, label='original')
plt.plot(t, s1(t),'o', label='knots')
plt.plot(x2, y2, '--', label='spline 1')
plt.plot(x2, s2(x2), 'r:', label='spline 2') 
plt.legend()

​

​

Answer 1

get_knots说：

在内部，节点向量包含2*k附加边界节点。

这意味着，要从get_knots获得一个可用的纽结数组，应该在数组的开头添加左边界结的k副本，在末尾添加右边界节点的k副本。这里，k是样条的度，通常是3(您要求LSQUnivariateSpline的缺省度，所以这是3)。所以：

kn = s1.get_knots()
kn = 3*[kn[0]] + list(kn) + 3*[kn[-1]]
c = s1.get_coeffs()
s2 = BSpline(kn, c, 3)    # not "2" as in your sample; we are working with a cubic spline 

现在，样条s2与s1相同：

​

​

同样地，kn = 4*[x[0]] + t + 4*[x[-1]]也能工作:您的t列表只包含内部结，因此添加了x[0]和x[-1]，然后每个重复k的次数都会更多。

重复的数学原因是B样条需要一些空间来建立，因为它们的归纳定义要求(k-1)-degree样条在定义k第四次样条的每一个区间附近存在。

Answer 2

这是一个稍微紧凑的方法，如果你不太关心结的细节位置。tk数组是您要寻找的。一旦tk在手中，样条可以使用y=splev(x,tk,der=0)线再现。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import splrep,splev
import matplotlib.pyplot as plt 

### Input data
x_arr = np.linspace(0, 10, 50)  # x-data
y_arr = np.exp(-(x_arr-5)**2/4)    # y-data
### Order of spline 
order = 3
### Make the spline 
tk = splrep(x_arr, y_arr, k=order) # Returns the knots and coefficents
### Evaluate the spline using the knots and coefficents on the domian x
x = np.linspace(0, 10, 1000) # new x-grid
y = splev(x, tk, der=0)
### Plot
f,ax=plt.subplots()
ax.scatter(x_arr, y_arr, label='original')
ax.plot(x,y,label='Spline')
ax.legend(fontsize=15)
plt.show()