我有几个如下所示的图:
我想知道有哪些方法可以找到 x 轴大约 5.5 到 8 之间的斜率。如果有几个这样的图,我更想知道是否有一种方法可以自动找到斜率值。
有什么建议么?
我在想ployfit()
,还是线性回归。问题是我不确定如何自动找到这些值。
我有几个如下所示的图:
我想知道有哪些方法可以找到 x 轴大约 5.5 到 8 之间的斜率。如果有几个这样的图,我更想知道是否有一种方法可以自动找到斜率值。
有什么建议么?
我在想ployfit()
,还是线性回归。问题是我不确定如何自动找到这些值。
在数据集中找到线性部分的一种通用方法是计算函数的二阶导数,并查看它在哪里(接近)为零。在解决问题的过程中需要考虑几件事:
如何计算噪声数据的二阶导数?一种快速简单的方法,可以很容易地适应不同的噪声水平、数据集大小和线性补丁的预期长度,是将数据与等于高斯二阶导数的卷积核进行卷积。可调整的部分是内核的宽度。
在您的上下文中,“接近于零”是什么意思?要回答这个问题,您必须对数据进行试验。
该方法的结果可以用作上述 chi^2 方法的输入,以识别数据集中的候选区域。
这里有一些可以帮助您入门的源代码:
from matplotlib import pyplot as plt
import numpy as np
# create theoretical data
x_a = np.linspace(-8,0, 60)
y_a = np.sin(x_a)
x_b = np.linspace(0,4,30)[1:]
y_b = x_b[:]
x_c = np.linspace(4,6,15)[1:]
y_c = np.sin((x_c - 4)/4*np.pi)/np.pi*4. + 4
x_d = np.linspace(6,14,120)[1:]
y_d = np.zeros(len(x_d)) + 4 + (4/np.pi)
x = np.concatenate((x_a, x_b, x_c, x_d))
y = np.concatenate((y_a, y_b, y_c, y_d))
# make noisy data from theoretical data
y_n = y + np.random.normal(0, 0.27, len(x))
# create convolution kernel for calculating
# the smoothed second order derivative
smooth_width = 59
x1 = np.linspace(-3,3,smooth_width)
norm = np.sum(np.exp(-x1**2)) * (x1[1]-x1[0]) # ad hoc normalization
y1 = (4*x1**2 - 2) * np.exp(-x1**2) / smooth_width *8#norm*(x1[1]-x1[0])
# calculate second order deriv.
y_conv = np.convolve(y_n, y1, mode="same")
# plot data
plt.plot(x,y_conv, label = "second deriv")
plt.plot(x, y_n,"o", label = "noisy data")
plt.plot(x, y, label="theory")
plt.plot(x, x, "0.3", label = "linear data")
plt.hlines([0],-10, 20)
plt.axvspan(0,4, color="y", alpha=0.2)
plt.axvspan(6,14, color="y", alpha=0.2)
plt.axhspan(-1,1, color="b", alpha=0.2)
plt.vlines([0, 4, 6],-10, 10)
plt.xlim(-2.5,12)
plt.ylim(-2.5,6)
plt.legend(loc=0)
plt.show()
这是结果:
smooth_width
是卷积核的宽度。为了调整噪声量,将0.27
random.normal 中的值更改为不同的值。请注意,此方法在靠近数据空间边界时效果不佳。
如您所见,二阶导数(蓝线)的“接近零”要求对于数据是线性的黄色部分非常适用。
您可以使用Ramer Douglas Peucker 算法将您的数据简化为一组较小的线段。该算法允许您指定一个epsilon
这样的方法,即每个数据点都不会比epsilon
某个线段更远。线段的斜率将粗略估计曲线的斜率。
这只是一种可能的解决方案,它会找到最小 chi^2 值大于预设最小值的点的直线段;
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()
我可以看到它对于更大的数据集会出现问题......
如果您的数据“模型”由大部分符合直线的数据组成,最后有一些异常值或摆动位,您可以尝试RANSAC算法。
这里的(非常罗嗦,抱歉)伪代码是:
choose a small threshold distance D
for N iterations:
pick two random points from your data, a and b
fit a straight line, L, to a and b
count the inliers: data points within a distance D of the line L
save the parameters of the line with the most inliers so far
estimate the final line using ALL the inliers of the best line