分段线性回归的优化

7

我试图创建一个分段线性回归来最小化均方差(MSE),然后直接使用线性回归。这个方法应该使用动态规划来计算不同分段大小和组合,以实现整体MSE的最小化。我认为算法运行时间复杂度为O(n²),我想知道是否有优化到O(nLogN)的方法?

import numpy as np
from sklearn.metrics import mean_squared_error
from sklearn import linear_model
import pandas as pd
import matplotlib.pyplot as plt

x = [3.4, 1.8, 4.6, 2.3, 3.1, 5.5, 0.7, 3.0, 2.6, 4.3, 2.1, 1.1, 6.1, 4.8,3.8]
y = [26.2, 17.8, 31.3, 23.1, 27.5, 36.5, 14.1, 22.3, 19.6, 31.3, 24.0, 17.3, 43.2, 36.4, 26.1]
dataset = np.dstack((x,y))
dataset = dataset[0]
d_arg = np.argsort(dataset[:,0])
dataset = dataset[d_arg]

def calc_error(dataset):
    lr_model = linear_model.LinearRegression()
    x = pd.DataFrame(dataset[:,0])
    y = pd.DataFrame(dataset[:,1])
    lr_model.fit(x,y)
    predictions = lr_model.predict(x)
    mse = mean_squared_error(y, predictions)
    return mse

#n is the number of points , m is the number of groups, k is the minimum number of points in a group
#(15,5,3)returns 【【3,3,3,3,3】】
#(15,5,2) returns [[2,4,3,3,3],[3,2,4,2,4],[4,2,3,3,3]....]
def all_combination(n,m,k):
    result = []
    if n < k*m:
        print('There are not enough elements to split.')
        return
    combination_bottom = [k for q in range(m)] 
    #add greedy algorithm here?
    if n == k*m:
        result.append(combination_bottom.copy())
    else:
        combination_now = [combination_bottom.copy()]
        j = k*m+1
        while j < n+1:
            combination_last = combination_now.copy()
            combination_now = []
            for x in combination_last:
                for i in range (0, m):
                    combination_new = x.copy()
                    combination_new[i] = combination_new[i]+1
                    combination_now.append(combination_new.copy())
            j += 1
        else:
            for x in combination_last:
                for i in range (0, m):
                    combination_new = x.copy()
                    combination_new[i] = combination_new[i]+1
                    if combination_new not in result:
                        result.append(combination_new.copy())
            
    return result #2-d list

def calc_sum_error(dataset,cb):#cb = combination
    mse_sum = 0    
    for n in range(0,len(cb)):
        if n == 0:
            low = 0
            high = cb[0]
        else:
            low = 0
            for i in range(0,n):
                low += cb[i]
            high = low + cb[n]
        mse_sum += calc_error(dataset[low:high])
    return mse_sum
            
    
    
#k is the number of points as a group
def best_piecewise(dataset,k):
    lenth = len(dataset)
    max_split = lenth // k
    min_mse = calc_error(dataset)
    split_cb = []
    all_cb = []
    for i in range(2, max_split+1):
        split_result = all_combination(lenth, i, k)
        all_cb += split_result
        for cb in split_result:
            tmp_mse = calc_sum_error(dataset,cb)
            if tmp_mse < min_mse:
                min_mse = tmp_mse
                split_cb = cb
    return min_mse, split_cb, all_cb

min_mse, split_cb, all_cb = best_piecewise(dataset, 2)

print('The best split of the data is '+str(split_cb))
print('The minimum MSE value is '+str(min_mse))

x = np.array(dataset[:,0])
y = np.array(dataset[:,1])

plt.plot(x,y,"o")
for n in range(0,len(split_cb)):
    if n == 0:
        low = 0 
        high = split_cb[n]
    else:
        low = 0
        for i in range(0,n):
            low += split_cb[i]
        high = low + split_cb[n]
    x_tmp = pd.DataFrame(dataset[low:high,0])
    y_tmp = pd.DataFrame(dataset[low:high,1])
    lr_model = linear_model.LinearRegression()
    lr_model.fit(x_tmp,y_tmp)
    y_predict = lr_model.predict(x_tmp)
    plt.plot(x_tmp, y_predict, 'g-')

plt.show()

如果有任何部分没有表述清楚,请告诉我。


你可能会发现Cubist回归模型很有用。该软件包提供了R绑定:https://cran.r-project.org/web/packages/Cubist/vignettes/cubist.html,您可以使用rpy2从Python中调用它。 - VersBersch
2个回答

7
我花了一些时间才意识到,你所描述的问题正是决策树回归器试图解决的问题。 不幸的是,构建一个最优的决策树是NP难问题,这意味着即使使用动态规划,你也无法将运行时间降至类似于O(NlogN)。
好消息是,你可以直接使用任何维护良好的决策树实现,例如sklearn.tree模块的DecisionTreeRegressor,并且可以确信在O(NlogN)时间复杂度下获得最佳性能。使用min_samples_leaf参数可以强制执行每个组的最小数据点数量。你还可以控制其他多种属性,如使用max_leaf_nodes控制最大组数,使用criterion优化不同损失函数等。
如果你想知道Scikit-learn的决策树与你的算法学到的决策树(即你代码中的split_cb)相比如何:
X = np.array(x).reshape(-1,1)
dt = DecisionTreeRegressor(min_samples_leaf=MIN_SIZE).fit(X,y)
split_cb = np.unique(dt.apply(X),return_counts=True)[1]

然后使用您通常使用的绘图代码。请注意,由于您的时间复杂度远高于O(NlogN)*,因此您的实现通常会找到比scikit-learn的贪心算法更好的拆分。


[1] Hyafil,L.,& Rivest,R. L.(1976)。构建最优二叉决策树是NP完全问题。信息处理信件,5(1),15-17

*尽管我不确定您的实现的确切时间复杂度,但肯定比O(N ^ 2)差,all_combination(21,4,2)耗时超过5分钟。


4

提醒一下,这是一个庞大的话题,无法在此讨论所有内容。但我认为我们可以通过回答你的问题来取得一些进展。

另外,我认为先了解理论最好,因为其他人可能没有到同样的地步。你的代码可以直接使用 - 太棒了! - 这有点表明你知道我要说什么,但还留下了一些问题:

  1. 为什么在不需要的情况下写纯Python,而且比NumPy慢得多,而你已经在某种程度上导入和使用它?

  2. 你的示例是否表明你不完全理解分段回归的应用?由于我们首先从理论开始,这可能不是一个问题。

关于回归的事情是这样的:它很少能完全模拟数据,越接近完美准确,就越容易过拟合。

你的分段回归,除了第一个之外,都是完美的。而且它们应该是这样的。两个点可以确定一条直线。因此,在你提供的示例中,你还展示了过度拟合数据和脆弱模型的例子。不确定是否正确吗?想想4.85到5.99的x值会返回什么?3.11到3.39呢?
你的示例在左侧(或顶部),标准的线性回归在右侧(或底部): 这个右侧的线性回归为我们提供了完整x值范围内的y值,并且(可能)会继续下去。函数的连续性正是我们所追求的。在另一个例子中,您可以尝试使用任意数量的工具,包括决策树回归器,但您要么得到类似脆弱的结果,要么得到不符合期望的结果。然后呢?因为它“错误”而抛弃它?因为“计算机说什么就是什么”而接受它?两种都一样糟糕。
我们可以在此停止。但这是一个非常好的问题,如果不继续探讨,那将是一种不公正。因此,让我们从两个不同的数据集开始,这些数据集我们知道是分段回归的好候选。
iterations = 500
x = np.random.normal(0, 1, iterations) * 10
y = np.where(x < 0, -4 * x + 3 , np.where(x < 10, x + 48, x + 98)) + np.random.normal(0, 3, iterations)

plt.scatter(x, y, s = 3, color = 'k')
plt.show()

...这给我们提供了左边的图像。我选择它有两个原因:它在x值上是连续的,但在y值上不是。右边的图像来自一个非常出色的R包,它在x轴上也是连续的,有一个明显的断点,但仍然最好通过三个分段回归来完成。稍后我会多提一些关于它的信息。

需要注意的几件事情:

  • 一种明显的检测断点的方法是查找哪里的线会通过单侧极限测试但不会通过双侧极限测试。也就是说,它不是可微的。

  • 在这样的虚拟数据中,检测断点会更加容易,因为我们使用代码来开发它们。但在现实生活中,我们可能需要一个不同的解决方案。让我们暂时搁置这个问题。

我之前提到的一个问题是想知道为什么你要编写原始的Python代码,而其他专门针对你的问题的库要快得多。那么,让我们找出有多快以及你可能会得到什么样的答案。并且让我们使用不连续的折磨测试来进行衡量。

from scipy import optimize

def piecewise_linear(x, x0, x1, b, k1, k2, k3):
    condlist = [x < x0, (x >= x0) & (x < x1), x >= x1]
    funclist = [lambda x: k1*x + b, lambda x: k1*x + b + k2*(x-x0), lambda x: k1*x + b + k2*(x-x0) + k3*(x - x1)]
    return np.piecewise(x, condlist, funclist)

p, e = optimize.curve_fit(piecewise_linear, x, y)

xd = np.linspace(-30, 30, iterations)
plt.plot(x, y, "ko" )
plt.plot(xd, piecewise_linear(xd, *p))

即使在这样一个相当极端的情况下,我们得到了一个快速而强大的答案,可能不太美观,需要考虑它是否是最佳的。因此,请花些时间来考虑这张图表。它是否是最优的(为什么)或不是最优的(为什么不是)?

enter image description here

顺便提一下,让我们来谈谈时间。在自己编写的版本上(导入,数据,绘图 - 整个过程)运行%%timeit的结果是:

10.8秒±160毫秒每次循环(7次运行的平均值±标准偏差,1次循环)

这比使用内置的NumPy和SciPy函数进行类似操作(但额外随机化500个数据点)慢了650倍。

16.5毫秒±1.87毫秒每次循环(7次运行的平均值±标准偏差,100次循环)

如果这对您来说还不够,这是非常合理的情况,因为(我有点透露底牌)我们期望分段线性回归可以捕捉到所有不连续的断点。因此,让我引荐您阅读datadog的GitHub Gist,因为a:没有必要重新发明轮子;b:他们有一个有趣的实现。随着代码的附带,还有一篇相关博客文章,介绍了动态规划的关键缺点以及他们的方法和思路。

虽然动态规划可用于比朴素的蛮力实现更有效地遍历此搜索空间,但在实践中仍然过慢。

最后三点。

  1. 如果你调整随机点以获得连续的线条,结果看起来更好,但仍不完美。
  2. 你可以看到为什么这不能在一个问题中全部解决。我没有涉及到像曲线拟合、安斯康姆四重奏、样条或使用 ML 等内容。
  3. 目前,没有什么能替代理解正在发生的事情。话虽如此,R 包 MCP 在使用贝叶斯方法识别拐点方面确实非常出色。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接