如何计算数据点在交叉点之外继续的方式？

Question

假设您有两个来自计算的数据值数组,您可以使用连续,可微分函数进行建模.数据点的两条"线"在(至少)一点交叉,现在的问题是这些数据集背后的函数是否实际交叉或反交叉.

下面的图片显示了这种情况,我知道(从它背后的物理学),在上面的两个"接触点",黄色和绿色线实际上应该"切换颜色",而在较低的一个,两个功能彼此不同办法:

为了给出更简单的"玩具集"数据,请以此代码为例:

import matplotlib.pyplot as plt
import numpy as np

x=np.arange(-10,10,.5)
y1=[np.absolute(i**3)+100*np.absolute(i) for i in x]
y2=[-np.absolute(i**3)-100*np.absolute(i) for i in x][::-1]

plt.scatter(x,y1)
plt.scatter(x,y2,color='r')

plt.show()

哪个应该产生以下图像:

现在我怎么能推断数据背后的趋势是否交叉(使得从数据下部左侧延续到上右)或抗交叉(如与上面的颜色表示,从数据下部左侧延续到下右)？

到目前为止,我能够通过查看它们之间差异的导数找到这些数据集之间的"联系点",大致如下:

closePoints=np.where(np.diff(np.diff(array_A - array_B) > 0))[0] + 1 

(用scipy的cKDTree来评估可能会更快).

我应该继续(可能非常低效)检查交叉路口两侧的导数吗？或者我可以以某种方式检查左侧数据的推断是否更适合穿越或反交叉？

Answer 1

您可以对差异函数使用样条插值g(x) = y1(x) - y(2)。找到平方的最小值g(x)**2将是接触点或交叉点。查看一阶和二阶导数，您可以确定它是接触点（g(x)具有最小值，g'(x)==0, g''(x) != 0）还是交叉点（g(x)是静止点，g'(x)==0, g''(x)==0）。

以下代码搜索约束区间内的最小值g(x)**2，然后绘制导数。约束区间的使用是通过排除先前点所在的区间来连续找到多个点。

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as sopt
import scipy.interpolate as sip

# test functions:
nocrossingTest = True
if nocrossingTest:
    f1 = lambda x: +np.absolute(x**3)+100*np.absolute(x)
    f2 = lambda x: -np.absolute(x**3)-100*np.absolute(x)
else:
    f1 = lambda x: +np.absolute(x**3)+100*x
    f2 = lambda x: -np.absolute(x**3)-100*x

xp = np.arange(-10,10,.5)
y1p, y2p = f1(xp), f2(xp) # test array


# Do Interpolation of y1-y2 to find crossing point:
g12 = sip.InterpolatedUnivariateSpline(xp, y1p - y2p) # Spline Interpolator of Difference
dg12 = g12.derivative() # spline derivative
ddg12 = dg12.derivative() # spline derivative

# Bounded least square fit to find minimal distance
gg = lambda x: g12(x)*g12(x)
rr = sopt.minimize_scalar(gg, bounds=[-1,1]) # search minium in Interval [-1,1]
x_c = rr['x'] # x value with minimum distance
print("Crossing point is at x = {} (Distance: {})".format(x_c, g12(x_c)))



fg = plt.figure(1)
fg.clf()
fg,ax = plt.subplots(1, 1,num=1)
ax.set_title("Function Values $y$")
ax.plot(xp, np.vstack([y1p,y2p]).T, 'x',)
xx = np.linspace(xp[0], xp[-1], 1000)
ax.plot(xx, np.vstack([f1(xx), f2(xx)]).T, '-', alpha=0.5)
ax.grid(True)
ax.legend(loc="best")
fg.canvas.draw()

fg = plt.figure(2)
fg.clf()
fg,axx = plt.subplots(3, 1,num=2)
axx[0].set_title("$g(x) = y_1(x) - y_2(x)$")
axx[1].set_title("$dg(x)/dx$")
axx[2].set_title("$d^2g(x)/dx^2$")
for ax,g in zip(axx, [g12, dg12, ddg12]):
    ax.plot(xx, g(xx))
    ax.plot(x_c, g(x_c), 'ro', alpha=.5)
    ax.grid(True)

fg.tight_layout()
plt.show()

差分函数表明差分不平滑：