使用 scipy 查找样条曲线的平滑度
Finding the smoothness of a spline using scipy
考虑以下示例:
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy import interpolate
xs = np.linspace(1,10,500)
ys = [0.92 * x ** 2.3 + 0.0132 * x ** 4 + 0.0743 * (x - 9) ** 3 - 4 * (x -3) ** 2 + 80 * math.sin(math.sin(x)) + 10 * math.sin(x*5) + 1.2* np.random.normal(-4,4,1) for x in xs]
ys[200] = ys[200] + 130
ys[201] = ys[201] + 135
ys[202] = ys[202] + 129
ys[203] = ys[203] + 128
ys[204] = ys[204] + 131
ys[205] = ys[205] + 130
ys[206] = ys[206] + 129
ys[207] = ys[207] + 129
ys[208] = ys[208] + 128
ys[209] = ys[209] + 130
如果我在此时绘制 xs
和 ys
,它会生成一个漂亮的图形:
现在我正在使用 scipy.interpolate.splrep
将样条曲线拟合到该数据。我使用了两个不同的样条曲线来拟合两个不同的数据段:
tck = interpolate.splrep(xs[0:199], ys[0:199], s = 1000)
ynew2 = interpolate.splev(xs[0:199], tck, der = 0)
和:
tck = interpolate.splrep(xs[210:500], ys[210:500], s = 9000)
ynew3 = interpolate.splev(xs[210:500], tck, der = 0)
然后我们有:
现在我想以编程方式检测拟合质量。拟合既不应该太直 - 即保留特征,也不应该 "overdetect" 噪声变化作为特征。
我计划使用馈送到 ANN 的峰值计数器。
但是,此时,我的问题是:
- scipy/numpy 是否有一个内置函数,我可以在其中输入
splrep
的输出,它会在任何特定时间间隔计算 maxima/minima 的最小值或最大值和密度?
注:
我知道 R**2
值,我正在寻找另一种方法来检测特征的保留。
SciPy 没有找到三次样条曲线的临界点的方法。我们有最接近的 sproot 可以找到三次样条的根。为了让它在这里有用,我们必须拟合 4 阶样条,以便导数是三次样条。这就是我在下面所做的
from scipy.interpolate import splrep, splev, splder, sproot
tck1 = splrep(xs[0:199], ys[0:199], k=4, s=1000)
tck2 = splrep(xs[210:500], ys[210:500], k=4, s=9000)
roots1 = sproot(splder(tck1), 1000) # 1000 is an upper bound for the number of roots
roots2 = sproot(splder(tck2), 1000)
x1 = np.linspace(xs[0], xs[198], 1000) # plot both splines
plt.plot(x1, splev(x1, tck1))
x2 = np.linspace(xs[210], xs[499], 1000)
plt.plot(x2, splev(x2, tck2))
plt.plot(roots1, splev(roots1, tck1), 'ro') # plot their max/min points
plt.plot(roots2, splev(roots2, tck2), 'ro')
plt.show()
区别很明显。
您还可以找到任何特定区间的根数,例如 [3, 4]:
np.where((3 <= roots1) & (roots1 <= 4))[0].size # 29
或等效地,np.sum((3 <= roots1) & (roots1 <= 4))
考虑以下示例:
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy import interpolate
xs = np.linspace(1,10,500)
ys = [0.92 * x ** 2.3 + 0.0132 * x ** 4 + 0.0743 * (x - 9) ** 3 - 4 * (x -3) ** 2 + 80 * math.sin(math.sin(x)) + 10 * math.sin(x*5) + 1.2* np.random.normal(-4,4,1) for x in xs]
ys[200] = ys[200] + 130
ys[201] = ys[201] + 135
ys[202] = ys[202] + 129
ys[203] = ys[203] + 128
ys[204] = ys[204] + 131
ys[205] = ys[205] + 130
ys[206] = ys[206] + 129
ys[207] = ys[207] + 129
ys[208] = ys[208] + 128
ys[209] = ys[209] + 130
如果我在此时绘制 xs
和 ys
,它会生成一个漂亮的图形:
现在我正在使用 scipy.interpolate.splrep
将样条曲线拟合到该数据。我使用了两个不同的样条曲线来拟合两个不同的数据段:
tck = interpolate.splrep(xs[0:199], ys[0:199], s = 1000)
ynew2 = interpolate.splev(xs[0:199], tck, der = 0)
和:
tck = interpolate.splrep(xs[210:500], ys[210:500], s = 9000)
ynew3 = interpolate.splev(xs[210:500], tck, der = 0)
然后我们有:
现在我想以编程方式检测拟合质量。拟合既不应该太直 - 即保留特征,也不应该 "overdetect" 噪声变化作为特征。
我计划使用馈送到 ANN 的峰值计数器。
但是,此时,我的问题是:
- scipy/numpy 是否有一个内置函数,我可以在其中输入
splrep
的输出,它会在任何特定时间间隔计算 maxima/minima 的最小值或最大值和密度?
注:
我知道 R**2
值,我正在寻找另一种方法来检测特征的保留。
SciPy 没有找到三次样条曲线的临界点的方法。我们有最接近的 sproot 可以找到三次样条的根。为了让它在这里有用,我们必须拟合 4 阶样条,以便导数是三次样条。这就是我在下面所做的
from scipy.interpolate import splrep, splev, splder, sproot
tck1 = splrep(xs[0:199], ys[0:199], k=4, s=1000)
tck2 = splrep(xs[210:500], ys[210:500], k=4, s=9000)
roots1 = sproot(splder(tck1), 1000) # 1000 is an upper bound for the number of roots
roots2 = sproot(splder(tck2), 1000)
x1 = np.linspace(xs[0], xs[198], 1000) # plot both splines
plt.plot(x1, splev(x1, tck1))
x2 = np.linspace(xs[210], xs[499], 1000)
plt.plot(x2, splev(x2, tck2))
plt.plot(roots1, splev(roots1, tck1), 'ro') # plot their max/min points
plt.plot(roots2, splev(roots2, tck2), 'ro')
plt.show()
区别很明显。
您还可以找到任何特定区间的根数,例如 [3, 4]:
np.where((3 <= roots1) & (roots1 <= 4))[0].size # 29
或等效地,np.sum((3 <= roots1) & (roots1 <= 4))