<p>在数据集中寻找线性部分的一种通用方法是计算函数的二阶导数,并查看它在哪里(接近)为零。在寻求解决方案的过程中,需要考虑以下几点:</p>
<ul>
<li><p>如何计算噪声数据的二阶导数?一种快速而简单的方法是用一个等于高斯二阶导数的卷积核对数据进行卷积,这种方法可以很容易地适应不同的噪声水平、数据集大小和线性面片的预期长度。可调部分是内核的宽度。</p></li>
<li><p>在你的上下文中,“接近零”是什么意思?要回答这个问题,你必须用你的数据做实验。</p></li>
<li><p>该方法的结果可以作为上述chi^2方法的输入,以识别数据集中的候选区域。</p></li>
</ul>
<p>下面是一些可以帮助您开始的源代码:</p>
<pre><code>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()
</code></pre>
<p>结果是:
<img src="https://i.stack.imgur.com/9KVs6.png" alt="enter image description here"/></p>
<p><code>smooth_width</code>是卷积核的宽度。为了调整噪声量,将random.normal中的值<code>0.27</code>更改为不同的值。请注意,这种方法在靠近数据空间边界的地方不太有效。</p>
<p>如您所见,二阶导数(蓝线)的“接近零”要求对黄色部分适用得很好,其中数据是线性的。</p>