假设我有
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
,使其重现样条曲线拟合的形状。我找到了here帖子,但是对于上面的示例,在那里产生的样条曲线似乎不正确,并且它也不正是我所需要的,因为它不返回表达式。
如何将样条曲线变成分段函数?是否有(简单)方式来表达每个间隔,例如作为多项式?
最佳答案
简短的答案是,如果您对基于标准幂的多项式系数感兴趣,那么最好使用 CubicSpline
(请参见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])
LONGE ANSWER:
这里有一些潜在的困惑点:
UnivariateSpline
适合B-spline基础,因此系数与标准多项式幂基础PPoly.from_spline
,但是不幸的是,UnivariateSpline
返回一个截断的结和系数列表,该函数无法使用。我们可以通过访问样条线对象的内部数据来解决此问题,这有点禁忌。 c
(无论来自UnivariateSpline
还是CubicSpline
)的阶次相反,并假设您以自己为“中心”,例如c[k,i]
处的系数属于c[k,i]*(x-x[i])^(3-k)
。 给定您的设置,请注意,如果不使用UnivariateSpline包装器,而直接使用
splrep
且不进行平滑处理(s=0
),则可以获取tck
(knots-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事情(运行所有这些的算法)。如果尝试从UnivariateSpline发送带有
tck
的(s.get_knots(), s.get_coeffs(), 3)
元组,则会丢失这些重复,因此from_spline
不起作用。 checkout 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
希望coeff的顺序相反,因此我们可以按原样传递p.c
。
关于python - 如何将样条曲线拟合转换为分段函数?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/61488797/