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

Question

我有一组数据点，看起来有点像一条直线，在开始处有一条曲线。请参见下图，其中显示了具有最佳拟合线的点（适合整个数据集）。在

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

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

我的一个想法是选择一个“cutouff point”，比如在t=100（x轴），然后将左边的点拟合成一条直线，将右边的点拟合成另一条线。然后我可以计算两条直线的chi^2，以得到拟合的“优度”。然后，我可以继续做同样的事情很多次，使用稍微不同的截止点，直到我找到一对线，给出最好的整体拟合。在

另一种更复杂的想法是将这些数据点描述为两条线的组合，y= m1*t + m2*t + b1 + b2，其中{}s是斜率，bs是y截距。然后，取总曲线的导数，我会得到dy/dt = m1+m2。然后也许我可以循环使用不同的“剪切点”，并拟合直线，直到得到一个导数最接近m1+m2的组合。但我不知道该怎么做，因为一开始我不是在处理一条曲线，而是一堆离散的数据点。在

优化Python中的两个fits的最佳方法是什么？在

Answer 1

下面是一个将两条直线拟合到一个数据集的示例，两条直线之间的交叉点也作为参数进行拟合。这个例子使用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)

Answer 2

这可以解释为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)))

然后我用^{}做两个线性拟合，它本身基于least squares的方法：

^{pr2}$

绘制结果如下图：

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

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

Answer 3

从形式上讲，您提到的第二种方法基本上是尝试将数据拟合到圆锥曲线上。本质上，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，它将非常适合您的数据集。在

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

这个方法在数学上非常简洁，但是我敢打赌，让你训练过的模型收敛到一个POSL，相当于在刀口上平衡某个东西（POSL的条件本质上是非常狭窄的）。最有可能的是，它最终会给你一个抛物线、椭圆或双曲线的方程（这可能只会使你的数据集得到充分的拟合，即使是非最优的二次曲线回归解也是值得的）。在

尽管如此，如果你不觉得结果令人满意，那么一个更好的方法可能就是不再担心形状，使用神经网络来进行这种形式的非线性回归。或者你可以观察肘关节点，然后按照你在第一种方法中建议的那样划分你的数据集。在