如何将两个线性函数拟合到一组数据点？

Question

我有一组数据点，它们看起来有点像一条在起点附近有一条曲线的线。请参见下图，它显示了具有最佳拟合线的点（适合整个数据集）。

相反，它们可以用两个线性函数来描述（一条穿过最左边的一组点的线和一条穿过其余数据点的单独的线）。这些点实际对应的是中子衰变，其中包含两种不同的同位素。

我不知道哪个点对应哪个同位素，所以我需要以某种方式做出最好的猜测。一种同位素的曲线将是一条直线，而另一种同位素的曲线将是一条不同的直线。如何将两条不同的最佳拟合线（线性）拟合到数据点集，从而优化两者的拟合？

我的一个想法是选择一个“截止点”，例如t=100（x 轴），并将左侧的点拟合为一条线，将右侧的点拟合为另一条线。然后我可以计算两条线的 chi^2 以获得拟合的“优点”。然后，我可以多次使用略有不同的截止点做同样的事情，直到找到能够提供最佳整体拟合的那对线。

另一个似乎更复杂的想法是将这些数据点描述为两条线的组合，y= m1*t + m2*t + b1 + b2，其中ms 是斜率，bs 是 y 截距。然后，取总曲线的导数，我会得到dy/dt = m1+m2。然后也许我可以循环通过不同的“截止点”，并拟合线，直到我得到一个导数最接近m1+m2 的组合。但我不知道该怎么做，因为我最初没有使用一条曲线，只是一堆离散的数据点。

在 Python 中优化两种拟合的最佳方法是什么？

Answer 1

这可以解释为time-series segmentation 与linear regression 结合的问题。有多种方法可以解决这个问题。您已经提到的其中之一：手动选择的数据分段点，另一个正在尝试最小化错误。

首先我尝试重新创建数据：

import numpy as np; import matplotlib.pyplot as plt
y1 = np.linspace(5.5, 3.7, num=100)
y1 = y1 + np.random.rand(y1.shape[0]) * np.linspace(0, .3, num=y1.shape[0])
y2 = np.linspace(3.7, 1.1, num=500)
y2 = y2 + np.random.rand(y2.shape[0]) * np.linspace(0.3, 1.9, num=y2.shape[0])
y = np.append(y1, y2)
x = np.array(range(len(y)))

然后我使用numpy.linalg.lstsq 进行两次线性拟合，这本身是基于least squares 的方法：

x1 = x[:100]
y1 = y[:100]
A1 = np.vstack([x1, np.ones(len(x1))]).T
m1, c1 = np.linalg.lstsq(A1, y1, rcond=None)[0]

x2 = x[100:]
y2 = y[100:]
A2 = np.vstack([x2, np.ones(len(x2))]).T
m2, c2 = np.linalg.lstsq(A2, y2, rcond=None)[0]

绘制结果如下图：

plt.scatter(x, y)
plt.plot(x1, m1 * x1 + c1, 'r')
plt.plot(x2, m2 * x2 + c2, 'r')
plt.show()

您现在可以使用 SWAB 之类的自动分割算法来替换 [100:] 和 [:100] 切片，但我建议您不要这样做，如果您能够手动决定何时拆分数据.如果您正在寻找实现，请查看提供的链接。

Answer 2

形式上，您提到的第二种方法本质上是尝试将您的数据拟合到圆锥曲线上。本质上，a pair of straight lines (POSL) is a type of conic section，因此可以用一般圆锥方程ax^2 + by^2 + 2hxy + 2gx + 2fy + c = 0 表示（而不是你在问题中提到的那个，它最终只是一条斜率 m1+m2 的直线）。为了澄清，肯定将有一个上述形式的方程，当绘制时，它会给出一个适合您的数据的 POSL。你的算法能否收敛到它是另一回事。

您可以尝试使用方法来找到系数 a、b、h、g、f 和 c。理想情况下，您为圆锥曲线获得的系数将形成一个 POSL，这将非常适合您的数据集。

如果您决定实现这一点，则必须记住，这个通用方程可以表示很多形状，如抛物线、双曲线等。训练后您可能会发现回归卡住了，或不收敛并采取不同的形状，如抛物线。您可以尝试推动回归到 POSL 形状，方法是奖励遵守these conditions 必要的圆锥截面成为回归方法中的 POSL。不过，这可能会使事情过于复杂。

这种方法在数学上非常简洁，但我敢打赌，试图让您训练的模型收敛到 POSL 将等同于在刀刃上平衡某些东西（POSL 的条件本质上非常狭窄）。最有可能的是，它最终会为您提供抛物线、椭圆或双曲线方程（这可能最终会充分拟合您的数据集，即使是非最优的圆锥回归解决方案也值得）。

不过，如果您对结果不满意，那么更好的方法可能是停止担心形状并使用神经网络进行这种形式的非线性回归。或者您可以观察肘点，然后按照您在第一种方法中的建议划分关于它的数据集。

Answer 3

这是将两条直线拟合到一个数据集的示例，两条直线之间的交叉点也作为参数拟合。此示例使用 scipy 的差分进化 (DE) 遗传算法来确定初始参数估计。 DE 的 scipy 实现使用拉丁超立方算法来确保对参数空间的彻底搜索，并且该算法需要搜索范围 - 在本示例中，这些范围取自数据最大值和最小值。

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings

xData = numpy.array([19.1647, 18.0189, 16.9550, 15.7683, 14.7044, 13.6269, 12.6040, 11.4309, 10.2987, 9.23465, 8.18440, 7.89789, 7.62498, 7.36571, 7.01106, 6.71094, 6.46548, 6.27436, 6.16543, 6.05569, 5.91904, 5.78247, 5.53661, 4.85425, 4.29468, 3.74888, 3.16206, 2.58882, 1.93371, 1.52426, 1.14211, 0.719035, 0.377708, 0.0226971, -0.223181, -0.537231, -0.878491, -1.27484, -1.45266, -1.57583, -1.61717])
yData = numpy.array([0.644557, 0.641059, 0.637555, 0.634059, 0.634135, 0.631825, 0.631899, 0.627209, 0.622516, 0.617818, 0.616103, 0.613736, 0.610175, 0.606613, 0.605445, 0.603676, 0.604887, 0.600127, 0.604909, 0.588207, 0.581056, 0.576292, 0.566761, 0.555472, 0.545367, 0.538842, 0.529336, 0.518635, 0.506747, 0.499018, 0.491885, 0.484754, 0.475230, 0.464514, 0.454387, 0.444861, 0.437128, 0.415076, 0.401363, 0.390034, 0.378698])


def func(xArray, breakpoint, slopeA, offsetA, slopeB, offsetB):
    returnArray = []
    for x in xArray:
        if x < breakpoint:
            returnArray.append(slopeA * x + offsetA)
        else:
            returnArray.append(slopeB * x + offsetB)
    return returnArray


# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = func(xData, *parameterTuple)
    return numpy.sum((yData - val) ** 2.0)


def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)
    slope = 10.0 * (maxY - minY) / (maxX - minX) # times 10 for safety margin

    parameterBounds = []
    parameterBounds.append([minX, maxX]) # search bounds for breakpoint
    parameterBounds.append([-slope, slope]) # search bounds for slopeA
    parameterBounds.append([minY, maxY]) # search bounds for offsetA
    parameterBounds.append([-slope, slope]) # search bounds for slopeB
    parameterBounds.append([minY, maxY]) # search bounds for offsetB


    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()

# call curve_fit without passing bounds from genetic algorithm
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Parameters:', fittedParameters)
print()

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))

print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    axes.plot(xData, yData,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(xData), max(xData))
    yModel = func(xModel, *fittedParameters)

    # now the model as a line plot
    axes.plot(xModel, yModel)

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)