为了解决您的问题，一个好的方法是定义一个与您的数据集的假设相匹配的概率模型。在您的情况下，您可能需要混合线性回归模型。您可以通过将不同的数据点与不同的混合分量相关联来创建类似于高斯混合模型的“回归量混合”模型。

我已经包含了一些代码来帮助您入门。该代码为两个回归量的混合实现了一个 EM 算法（它应该相对容易扩展到更大的混合）。对于随机数据集，该代码似乎相当健壮。但是，与线性回归不同，混合模型具有非凸目标，因此对于真实数据集，您可能需要使用不同的随机起点运行一些试验。

import numpy as np
import matplotlib.pyplot as plt 
import scipy.linalg as lin

#generate some random data
N=100
x=np.random.rand(N,2)
x[:,1]=1

w=np.random.rand(2,2)
y=np.zeros(N)

n=int(np.random.rand()*N)
y[:n]=np.dot(x[:n,:],w[0,:])+np.random.normal(size=n)*.01
y[n:]=np.dot(x[n:,:],w[1,:])+np.random.normal(size=N-n)*.01


rx=np.ones( (100,2) )
r=np.arange(0,1,.01)
rx[:,0]=r

#plot the random dataset
plt.plot(x[:,0],y,'.b')
plt.plot(r,np.dot(rx,w[0,:]),':k',linewidth=2)
plt.plot(r,np.dot(rx,w[1,:]),':k',linewidth=2)

# regularization parameter for the regression weights
lam=.01

def em():
    # mixture weights
    rpi=np.zeros( (2) )+.5

    # expected mixture weights for each data point
    pi=np.zeros( (len(x),2) )+.5

    #the regression weights
    w1=np.random.rand(2)
    w2=np.random.rand(2)

    #precision term for the probability of the data under the regression function 
    eta=100

    for _ in xrange(100):
        if 0:
            plt.plot(r,np.dot(rx,w1),'-r',alpha=.5)
            plt.plot(r,np.dot(rx,w2),'-g',alpha=.5)

        #compute lhood for each data point
        err1=y-np.dot(x,w1)
        err2=y-np.dot(x,w2)
        prbs=np.zeros( (len(y),2) )
        prbs[:,0]=-.5*eta*err1**2
        prbs[:,1]=-.5*eta*err2**2

        #compute expected mixture weights
        pi=np.tile(rpi,(len(x),1))*np.exp(prbs)
        pi/=np.tile(np.sum(pi,1),(2,1)).T

        #max with respect to the mixture probabilities
        rpi=np.sum(pi,0)
        rpi/=np.sum(rpi)

        #max with respect to the regression weights
        pi1x=np.tile(pi[:,0],(2,1)).T*x
        xp1=np.dot(pi1x.T,x)+np.eye(2)*lam/eta
        yp1=np.dot(pi1x.T,y)
        w1=lin.solve(xp1,yp1)

        pi2x=np.tile(pi[:,1],(2,1)).T*x
        xp2=np.dot(pi2x.T,x)+np.eye(2)*lam/eta
        yp2=np.dot(pi[:,1]*y,x)
        w2=lin.solve(xp2,yp2)

        #max wrt the precision term
        eta=np.sum(pi)/np.sum(-prbs/eta*pi)

        #objective function - unstable as the pi's become concentrated on a single component
        obj=np.sum(prbs*pi)-np.sum(pi[pi>1e-50]*np.log(pi[pi>1e-50]))+np.sum(pi*np.log(np.tile(rpi,(len(x),1))))+np.log(eta)*np.sum(pi)
        print obj,eta,rpi,w1,w2

        try:
            if np.isnan(obj): break
            if np.abs(obj-oldobj)<1e-2: break
        except:
            pass

        oldobj=obj

    return w1,w2


#run the em algorithm and plot the solution
rw1,rw2=em()
plt.plot(r,np.dot(rx,rw1),'-r')
plt.plot(r,np.dot(rx,rw2),'-g')

plt.show()

在此线程的其他地方，user1149913 提供了很好的建议（定义概率模型）和强大方法的代码（EM 估计）。 有两个问题需要解决：

1. 如何应对偏离概率模型的情况（这在 2011-2012 年的数据中非常明显，在坡度较小的点的起伏中也有些明显）。

2. 如何确定 EM 算法（或任何其他算法）的良好起始值。

要解决 #2，请考虑使用Hough 变换。这是一种特征检测算法，用于寻找特征的线性延伸，可以有效地计算为Radon 变换。

从概念上讲，霍夫变换描绘了一组线。平面中的一条线可以通过它的斜率来参数化，$x$，以及它的距离，$y$，来自固定的原点。这其中的一个点$x,y$坐标系由此指定一条线。原始图中的每个点都确定了穿过该点的线铅笔：这铅笔在霍夫变换中显示为曲线。当原始图中的特征沿着一条公共线或足够接近一条线时，它们在霍夫变换中产生的曲线集合往往具有对应于该公共线的公共交点。通过在霍夫变换中找到这些强度最大的点，我们可以读出原始问题的良好解决方案。

为了开始使用这些数据，我首先剪掉了辅助材料（轴、刻度线和标签），并且为了更好地测量右下角明显离群的点，我首先剪掉了底部的轴。（当这些东西没有被裁剪掉时，程序仍然运行良好，但它还检测到轴、帧、刻度的线性序列、标签的线性序列，甚至是零星地位于底部轴上的点！）

img = Import["http://i.stack.imgur.com/SkEm3.png"]
i = ColorNegate[Binarize[img]]
crop2 = ImageCrop[ImageCrop[i, {694, 531}, {Left, Bottom}], {565, 467}, {Right, Top}]

（这个和其余的代码都在Mathematica中。）

此图像中的每个点对应于霍夫变换中的窄范围曲线，在此处可见。它们是正弦波：

hough2 = Radon[crop2, Method -> "Hough"]  // ImageAdjust

这在视觉上体现了问题是线聚类问题的意义：霍夫变换将其简化为点聚类问题，我们可以应用任何我们喜欢的聚类方法。

在这种情况下，聚类非常清晰，只需对霍夫变换进行简单的后处理即可。为了识别变换中强度最大的位置，我增加了对比度并在大约 1% 的半径范围内模糊了变换：这与原始图像中绘图点的直径相当。

blur = ImageAdjust[Blur[ImageAdjust[hough2, {1, 0}], 8]]

阈值化结果将其缩小为两个微小的斑点，其质心合理地确定了强度最大的点：这些估计了拟合线。

comp = MorphologicalComponents[blur, 0.777]) // Colorize

（阈值$0.777$凭经验发现：它只产生两个区域，两个区域中较小的一个几乎尽可能小。）

图像的左侧对应于 0 度（水平）的方向，当我们从左到右看时，该角度线性增加到 180 度。插值，我计算出这两个斑点分别以 19 度和 57.1 度为中心。我们还可以从 blob 的垂直位置读取截距。此信息产生初始拟合：

width = ImageDimensions[blur][[1]];
slopes =  Module[{x, y, z}, ComponentMeasurements[comp, "Centroid"] /. 
          Rule[x_, {y_, z_}] :>  Round[((y - 1/2)/(width - 1))  180., 0.1]
  ]

{19., 57.1}

以类似的方式，可以计算对应于这些斜率的截距，给出这些拟合：

（红线对应上图中的小粉红点，蓝线对应较大的水色斑点。）

在很大程度上，这种方法自动处理了第一个问题：与线性的偏差会抹掉强度最大的点，但通常不会使它们发生太大变化。坦率地说，异常点会在整个霍夫变换中产生低级噪声，这些噪声会在后处理过程中消失。

在这一点上，可以将这些估计值作为 EM 算法或似然最小化器的起始值（在给定好的估计值的情况下，将很快收敛）。不过，更好的是使用稳健的回归估计器，例如迭代重新加权最小二乘法。它能够为每个点提供回归权重。低权重表示一个点不“属于”一条线。如果需要，利用这些权重将每个点分配给其正确的线。然后，对点进行分类后，您可以对两组点分别使用普通最小二乘法（或任何其他回归程序）。

我发现这个问题与另一个问题有关。我实际上对这类问题进行了学术研究。请检查我的答案“最小平方根”拟合？具有多个最小值的拟合方法以获取更多详细信息。

whuber 的基于霍夫变换的方法对于您给出的简单场景来说是一个非常好的解决方案。我处理了具有更复杂数据的场景，例如：

我和我的合著者将此称为“数据关联”问题。当您尝试解决它时，主要问题通常是组合问题，因为可能的数据组合呈指数级。

我们有一个出版物“用于数据关联问题的高斯过程的重叠混合”，其中我们使用迭代技术处理了 N 曲线的一般问题，得到了非常好的结果。您可以在论文中找到链接的 Matlab 代码。

[更新] OMGP 技术的 Python 实现可以在GPClust 库中找到。

我有另一篇论文，我们放宽了问题以获得凸优化问题，但尚未被接受发表。它特定于 2 条曲线，因此它可以完美地处理您的数据。如果您有兴趣，请告诉我。

user1149913 有一个很好的答案 (+1)，但在我看来，您的数据收集在 2011 年底崩溃了，所以您必须切断这部分数据，然后仍然以不同的随机数运行几次开始系数，看看你得到什么。

一种直接的方法是通过肉眼将您的数据分成两组，然后使用您习惯的任何线性模型技术。在 R 中，它将是lm函数。

或者用眼睛拟合两条线。在 R 中，您将使用它abline来执行此操作。

数据混乱，有异常值，最后分崩离析，但 by-eye 有两条相当明显的线，所以我不确定花哨的方法是否值得。