如何将样条拟合转换为分段函数？

Question

假设我有

import numpy as np
from scipy.interpolate import UnivariateSpline

# "true" data; I don't know this function
x = np.linspace(0, 100, 1000)
d = np.sin(x * 0.5) + 2 + np.cos(x * 0.1)

# sample data; that's what I actually measured
x_sample = x[::20]
d_sample = d[::20]

# fit spline
s = UnivariateSpline(x_sample, d_sample, k=3, s=0.005)

plt.plot(x, d)
plt.plot(x_sample, d_sample, 'o')
plt.plot(x, s(x))
plt.show()

我明白了

我现在想要的是所有橙色点之间的函数，比如

knots = s.get_knots()
f0 = <some expression> for knots[0] <= x < knots[1]
f1 = <some expression> for knots[1] <= x < knots[2]
...

因此，fi 的选择方式应能重现样条拟合的形状。

我找到了 post here，但是上面的示例生成的样条曲线似乎不正确，而且它也不完全是我需要的，因为它没有 return 表达式。

如何将样条曲线转换为分段函数？是否有一种（直接的）方式来表达每个间隔，例如作为多项式？

Answer 1

应该可以像您描述的那样使用分段函数，例如：

import numpy as np
from scipy.interpolate import UnivariateSpline

# "true" data; I don't know this function
x = np.linspace(0, 100, 1000)
d = np.sin(x * 0.5) + 2 + np.cos(x * 0.1)

# sample data; that's what I actually measured
x_sample = x[::20]
d_sample = d[::20]

# fit spline
s = UnivariateSpline(x_sample, d_sample, k=3, s=0.005)

plt.plot(x, d)
plt.plot(x_sample, d_sample, 'o')
plt.plot(x, s(x))

knots = s.get_knots()

conditions = [x < knots[0], (x >= knots[0]) * (x < knots[1]), (x >= knots[1]) * (x < knots[10]), x >= knots[10]]
# need one function for each condition
fns = [0, lambda x :-x, lambda x: 0.01*x**2, lambda x: 0.2*x**0.5]
y = np.piecewise(x, conditions, fns)

plt.plot(x, y)
plt.show()

我选择了一些随机的条件和函数-我相信你能找到更合适的！

Answer 2

简短的回答是，如果您对标准幂基多项式的系数感兴趣，最好使用 CubicSpline (and see this discussion):

cu = scipy.interpolate.CubicSpline(x_sample, d_sample)

plt.plot(x_sample, y_sample, 'ko')
for i in range(len(cu.x)-1):
    xs = np.linspace(cu.x[i], cu.x[i+1], 100)
    plt.plot(xs, np.polyval(cu.c[:,i], xs - cu.x[i]))

为了回答您的问题，您可以使用 numpy.piecewise、cu.x 中的断点和 cu.c 中的系数从此处创建一个分段函数，然后直接对多项式表达式自己或使用 numpy.polyval。例如，

cu.c[:,0]  # coeffs for 0th segment
# array([-0.01316353, -0.02680068,  0.51629024,  3.        ])

# equivalent ways to code polynomial for this segment
f0 = lambda x: cu.c[0,0]*(x-x[0])**3 + cu.c[1,0]*(x-x[0])**2 + cu.c[2,0]*(x-x[0]) + cu.c[3,0]
f0 = lambda x: [cu.c[i,0]*(x-x[0])**(3-i) for i in range(4)]

# ... or getting values directly from x's
y0 = np.polyval(cu.c[:,0], xs - cu.x[0])

更长的答案：

这里有几点可能引起混淆：

UnivariateSpline 适合 B-spline 基，因此系数与标准多项式幂基
为了从 B 样条转换，我们可以使用 PPoly.from_spline，但不幸的是 UnivariateSpline returns 一个 被截断的 节点列表和不会与此功能一起玩的系数。我们可以通过访问样条对象的内部数据来解决这个问题，这个有点忌讳。
此外，系数矩阵 c（无论是来自 UnivariateSpline 还是 CubicSpline）的度数倒序并假定您是 "centering" 自己，例如c[k,i] 处的系数属于 c[k,i]*(x-x[i])^(3-k).

根据您的设置，请注意，如果我们不使用 UnivariateSpline 包装器，而是直接使用 splrep 进行拟合并且没有平滑 (s=0)，我们可以获取 tck（节-coefficients-degree) 元组并将其发送到 PPoly.from_spline 函数并得到我们想要的系数：

tck = scipy.interpolate.splrep(x_sample, d_sample, s=0)
tck
# (array([0.        , 0.        , 0.        , 0.        , 2.68456376,
#        4.02684564, 5.36912752, 6.7114094 , 9.39597315, 9.39597315,
#        9.39597315, 9.39597315]),
# array([3.        , 3.46200469, 4.05843704, 3.89649312, 3.33792889,
#        2.29435138, 1.65015175, 1.59021688, 0.        , 0.        ,
#        0.        , 0.        ]),
# 3)

p = scipy.interpolate.PPoly.from_spline(tck)
p.x.shape  # breakpoints in unexpected shape
# (12,)

p.c.shape  # polynomial coeffs in unexpected shape
# (4, 11)

请注意 tck 和 p.x 中奇怪的重复断点：这是一个 FITPACK 的东西（算法运行所有这些）。

如果我们尝试使用 (s.get_knots(), s.get_coeffs(), 3) 从 UnivariateSpline 发送一个 tck 元组，我们会丢失那些重复项，因此 from_spline 不起作用。查看 source 虽然看起来完整的向量存储在 self._data 中，所以我们可以做

s = scipy.interpolate.UnivariateSpline(x_sample, d_sample, s=0)
tck = (s._data[8], s._data[9], 3)
p = scipy.interpolate.PPoly.from_spline(tck)

和以前一样。要检查这些系数是否有效：

plt.plot(x_sample, d_sample, 'o')

for i in range(len(p.x)-1):
    xs = np.linspace(p.x[i], p.x[i+1], 10)
    plt.plot(xs, np.polyval(p.c[:,i], xs - p.x[i]))

注意 numpy.polyval 需要 coeffs 的反向顺序，因此我们可以按原样传递 p.c。

如何将样条拟合转换为分段函数？

How to convert a spline fit into a piecewise function?

python

numpy

scipy

piecewise