from matplotlib.pyplot import figure, show
from numpy import pi, sin, linspace, exp, polyfit
from matplotlib.mlab import stineman_interp
x = linspace(0,2*pi,20);
y = x + sin(x) + exp(-0.5*(x-2)**2);
num_points = len(x)
min_fit_length = 5
chi = 0
chi_min = 10000
i_best = 0
j_best = 0
for i in range(len(x) - min_fit_length):
for j in range(i+min_fit_length, len(x)):
coefs = polyfit(x[i:j],y[i:j],1)
y_linear = x * coefs[0] + coefs[1]
chi = 0
for k in range(i,j):
chi += ( y_linear[k] - y[k])**2
if chi < chi_min:
i_best = i
j_best = j
chi_min = chi
print chi_min
coefs = polyfit(x[i_best:j_best],y[i_best:j_best],1)
y_linear = x[i_best:j_best] * coefs[0] + coefs[1]
fig = figure()
ax = fig.add_subplot(111)
ax.plot(x,y,'ro')
ax.plot(x[i_best:j_best],y_linear,'b-')
show()
在数据集中寻找线性部分的一种通用方法是计算函数的二阶导数,并查看它在哪里(接近)为零。在寻求解决方案的过程中,需要考虑以下几点:
如何计算噪声数据的二阶导数?一种快速而简单的方法是用一个等于高斯二阶导数的卷积核对数据进行卷积,这种方法可以很容易地适应不同的噪声水平、数据集大小和线性面片的预期长度。可调部分是内核的宽度。
在你的上下文中,“接近零”是什么意思?要回答这个问题,你必须用你的数据做实验。
该方法的结果可以作为上述chi^2方法的输入,以识别数据集中的候选区域。
下面是一些可以帮助您开始的源代码:
结果是:
smooth_width
是卷积核的宽度。为了调整噪声量,将random.normal中的值0.27
更改为不同的值。请注意,这种方法在靠近数据空间边界的地方不太有效。如您所见,二阶导数(蓝线)的“接近零”要求对黄色部分适用得很好,其中数据是线性的。
这只是一个可能的解决方案,它将找到点的直线段,其最小chi^2值大于预设的最小值
我可以看到,对于更大的数据集来说,这会变得有问题。。。
可以使用Ramer Douglas Peucker algorithm将数据简化为一组较小的线段。该算法允许您指定一个
epsilon
,这样每个数据点与某些线段的距离都不会超过epsilon
。线段的坡度可以粗略估计曲线的坡度。有一个Python implementation of the RDP algorithm here.
相关问题 更多 >
编程相关推荐