如何在 Python 中应用分段线性拟合？

Posted 2023-03-12

【中文标题】如何在 Python 中应用分段线性拟合？

【英文标题】：How to apply piecewise linear fit in Python?

【发布时间】：2015-06-05 15:13:33

【问题描述】：

我正在尝试对数据集进行分段线性拟合，如图 1 所示

这个数字是通过对行的设置得到的。我尝试使用代码应用分段线性拟合：

from scipy import optimize
import matplotlib.pyplot as plt
import numpy as np


x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15])
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03])


def linear_fit(x, a, b):
    return a * x + b
fit_a, fit_b = optimize.curve_fit(linear_fit, x[0:5], y[0:5])[0]
y_fit = fit_a * x[0:7] + fit_b
fit_a, fit_b = optimize.curve_fit(linear_fit, x[6:14], y[6:14])[0]
y_fit = np.append(y_fit, fit_a * x[6:14] + fit_b)


figure = plt.figure(figsize=(5.15, 5.15))
figure.clf()
plot = plt.subplot(111)
ax1 = plt.gca()
plot.plot(x, y, linestyle = '', linewidth = 0.25, markeredgecolor='none', marker = 'o', label = r'\textity_a')
plot.plot(x, y_fit, linestyle = ':', linewidth = 0.25, markeredgecolor='none', marker = '', label = r'\textity_b')
plot.set_ylabel('Y', labelpad = 6)
plot.set_xlabel('X', labelpad = 6)
figure.savefig('test.pdf', box_inches='tight')
plt.close()    

但这让我拟合了图 1 中的形式。 2，我尝试使用这些值，但没有改变我无法正确拟合上线。对我来说最重要的要求是如何让 Python 获得梯度变化点。本质上 我希望 Python 能够识别并拟合适当范围内的两个线性拟合。这如何在 Python 中完成？

【参考方案1】：

你可以使用numpy.piecewise()创建分段函数，然后使用curve_fit()，这里是代码

from scipy import optimize
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15], dtype=float)
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03])

def piecewise_linear(x, x0, y0, k1, k2):
    return np.piecewise(x, [x < x0], [lambda x:k1*x + y0-k1*x0, lambda x:k2*x + y0-k2*x0])

p , e = optimize.curve_fit(piecewise_linear, x, y)
xd = np.linspace(0, 15, 100)
plt.plot(x, y, "o")
plt.plot(xd, piecewise_linear(xd, *p))

输出：

N件配件请参考segments_fit.ipynb

【讨论】：

如何将其扩展到三个部分？

【参考方案2】：

您可以执行spline interpolation 方案来执行分段线性插值并找到曲线的转折点。二阶导数将在转折点处最高（对于单调递增的曲线），并且可以使用阶数 > 2 的样条插值来计算。

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15])
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03])

tck = interpolate.splrep(x, y, k=2, s=0)
xnew = np.linspace(0, 15)

fig, axes = plt.subplots(3)

axes[0].plot(x, y, 'x', label = 'data')
axes[0].plot(xnew, interpolate.splev(xnew, tck, der=0), label = 'Fit')
axes[1].plot(x, interpolate.splev(x, tck, der=1), label = '1st dev')
dev_2 = interpolate.splev(x, tck, der=2)
axes[2].plot(x, dev_2, label = '2st dev')

turning_point_mask = dev_2 == np.amax(dev_2)
axes[2].plot(x[turning_point_mask], dev_2[turning_point_mask],'rx',
             label = 'Turning point')
for ax in axes:
    ax.legend(loc = 'best')

plt.show()

【参考方案3】：

您可以使用 pwlf 在 Python 中执行连续分段线性回归。这个库可以使用 pip 安装。

pwlf 中有两种方法来执行您的匹配：

您可以适应指定数量的线段。 您可以指定连续分段线应终止的 x 位置。

让我们使用方法 1，因为它更容易，并且会识别您感兴趣的“梯度变化点”。

查看数据时，我注意到两个不同的区域。因此，使用两条线段找到可能的最佳连续分段线是有意义的。这是方法 1。

import numpy as np
import matplotlib.pyplot as plt
import pwlf

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15])
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59,
              84.47, 98.36, 112.25, 126.14, 140.03])

my_pwlf = pwlf.PiecewiseLinFit(x, y)
breaks = my_pwlf.fit(2)
print(breaks)

[1. 5.99819559 15.]

第一条线段从 [1., 5.99819559] 开始，第二条线段从 [5.99819559, 15.] 开始。因此，您要求的梯度变化点将是 5.99819559。

我们可以使用预测函数绘制这些结果。

x_hat = np.linspace(x.min(), x.max(), 100)
y_hat = my_pwlf.predict(x_hat)

plt.figure()
plt.plot(x, y, 'o')
plt.plot(x_hat, y_hat, '-')
plt.show()

【讨论】：

我刚刚发现了 pwlf，它似乎是一个非常好的插件包，有据可查

【参考方案4】：

两个变化点的示例。如果您愿意，只需根据此示例测试更多更改点即可。

np.random.seed(9999)
x = np.random.normal(0, 1, 1000) * 10
y = np.where(x < -15, -2 * x + 3 , np.where(x < 10, x + 48, -4 * x + 98)) + np.random.normal(0, 3, 1000)
plt.scatter(x, y, s = 5, color = u'b', marker = '.', label = 'scatter plt')

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, 1000)
plt.plot(x, y, "o")
plt.plot(xd, piecewise_linear(xd, *p))

【参考方案5】：

此方法使用Scikit-Learn 应用分段线性回归。 如果您的点受到噪音的影响，您可以使用它。 它比执行一项巨大的优化任务（来自scip.optimize 的任何东西，比如curve_fit 更多然后是 3 个参数）。

import numpy as np
import matplotlib.pylab as plt
from sklearn.tree import DecisionTreeRegressor
from sklearn.linear_model import LinearRegression

# parameters for setup
n_data = 20

# segmented linear regression parameters
n_seg = 3

np.random.seed(0)
fig, (ax0, ax1) = plt.subplots(1, 2)

# example 1
#xs = np.sort(np.random.rand(n_data))
#ys = np.random.rand(n_data) * .3 + np.tanh(5* (xs -.5))

# example 2
xs = np.linspace(-1, 1, 20)
ys = np.random.rand(n_data) * .3 + np.tanh(3*xs)

dys = np.gradient(ys, xs)

rgr = DecisionTreeRegressor(max_leaf_nodes=n_seg)
rgr.fit(xs.reshape(-1, 1), dys.reshape(-1, 1))
dys_dt = rgr.predict(xs.reshape(-1, 1)).flatten()

ys_sl = np.ones(len(xs)) * np.nan
for y in np.unique(dys_dt):
    msk = dys_dt == y
    lin_reg = LinearRegression()
    lin_reg.fit(xs[msk].reshape(-1, 1), ys[msk].reshape(-1, 1))
    ys_sl[msk] = lin_reg.predict(xs[msk].reshape(-1, 1)).flatten()
    ax0.plot([xs[msk][0], xs[msk][-1]],
             [ys_sl[msk][0], ys_sl[msk][-1]],
             color='r', zorder=1)

ax0.set_title('values')
ax0.scatter(xs, ys, label='data')
ax0.scatter(xs, ys_sl, s=3**2, label='seg lin reg', color='g', zorder=5)
ax0.legend()

ax1.set_title('slope')
ax1.scatter(xs, dys, label='data')
ax1.scatter(xs, dys_dt, label='DecisionTree', s=2**2)
ax1.legend()

plt.show()

它是如何工作的

计算每个点的斜率 使用决策树对相似的斜坡进行分组（右图） 对原始数据中的每一组进行线性回归

【讨论】：

这是一个非常好的方法。需要进行一些修改才能使其运行。 1：定义n_data和n_seg（并使用n_data生成xs）。 2：可能最好指定导入。

【参考方案6】：

使用numpy.interp 将一维分段线性插值返回给在离散数据点处具有给定值的函数。

【讨论】：

这个答案没有解决本质问题“我希望 Python 能够识别并拟合适当范围内的两个线性拟合。如何在 Python 中做到这一点？” numpy.interp 仅连接点，但不应用拟合。给出的示例的结果恰好是相同的，但通常情况并非如此。

【参考方案7】：

扩展@binoy-pilakkat 的答案。

你应该使用numpy.interp:

import numpy as np
import matplotlib.pyplot as plt

x = np.array(range(1,16), dtype=float)
y = np.array([5, 7, 9, 11, 13, 15, 28.92,
          42.81, 56.7, 70.59, 84.47,
          98.36, 112.25, 126.14, 140.03], dtype=float)

yinterp = np.interp(x, x, y) # simple as that

plt.plot(x, y, 'bo')
plt.plot(x, yinterp, 'g-')
plt.show()

【讨论】：

这个答案没有解决本质问题“我希望 Python 能够识别并拟合适当范围内的两个线性拟合。如何在 Python 中做到这一点？” numpy.interp 仅连接点，但不应用拟合。给出的示例的结果恰好是相同的，但通常情况并非如此。

【参考方案8】：

piecewise 也可以使用

from piecewise.regressor import piecewise
import numpy as np

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15,16,17,18], dtype=float)
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03,120,112,110])

model = piecewise(x, y)

评估“模型”：

FittedModel with segments:
* FittedSegment(start_t=1.0, end_t=7.0, coeffs=(2.9999999999999996, 2.0000000000000004))
* FittedSegment(start_t=7.0, end_t=16.0, coeffs=(-68.2972222222222, 13.888333333333332))
* FittedSegment(start_t=16.0, end_t=18.0, coeffs=(198.99999999999997, -5.000000000000001))

【参考方案9】：

我认为scipy.interpolate 中的UnivariateSpline 将提供最简单且很可能是最快的分段拟合方法。为了添加一些上下文，样条是由多项式分段定义的函数。在您的情况下，您正在寻找由k=1 在UnivariateSpline 中定义的线性样条曲线。此外，s=0.5 是一个平滑因子，表示拟合应该有多好（查看文档以获取更多信息）。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15])
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03])


# Solution
spl = UnivariateSpline(x, y, k=1, s=0.5)
xs = np.linspace(x.min(), x.max(), 1000)


fig, ax = plt.subplots()
ax.scatter(x, y, color="red", s=20, zorder=20)
ax.plot(xs, spl(xs), linestyle="--", linewidth=1, color="blue", zorder=10)
ax.grid(color="grey", linestyle="--", linewidth=.5, alpha=.5)
ax.set_ylabel("Y")
ax.set_xlabel("X")
plt.show()

【参考方案10】：

这里已经有很好的答案，但这里有另一种使用简单神经网络的方法。基本思想与其他一些答案相同；即，

创建指示输入变量是否大于某个断点的虚拟变量 通过从输入变量中减去断点然后将结果与相应的虚拟变量相乘来创建虚拟交互 使用输入变量和虚拟交互作为特征来训练线性模型

主要区别在于，这里的断点是通过梯度下降端到端学习的，而不是被视为超参数。这种方法自然会扩展到多个断点，并且可以与任何相关的损失函数一起使用。

import torch
import numpy as np
import matplotlib.pyplot as plt

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15])
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59,
              84.47, 98.36, 112.25, 126.14, 140.03])

定义模型、优化器和损失函数：

class PiecewiseLinearModel(torch.nn.Module):
    def __init__(self, n_breaks):
        super(PiecewiseLinearModel, self).__init__()
        self.breaks = torch.nn.Parameter(torch.randn((1,n_breaks)))
        self.linear = torch.nn.Linear(n_breaks+1, 1)        
    def forward(self, x):
        return self.linear(torch.cat([x, torch.nn.ReLU()(x - self.breaks)],1))

plm = PiecewiseLinearModel(n_breaks=1)
optimizer = torch.optim.Adam(plm.parameters(), lr=0.1)
loss_func = torch.nn.functional.mse_loss

训练模型：

x_torch = torch.tensor(x, dtype=torch.float)[:,None]
y_torch = torch.tensor(y)[:,None]
for _ in range(10000):
    p = plm(x_torch)
    optimizer.zero_grad()
    loss_func(y_torch, p).backward()
    optimizer.step()

绘制预测：

x_grid = np.linspace(0,16,1000)
p = plm(torch.tensor(x_grid, dtype=torch.float)[:,None])
p = p.flatten().detach().numpy()
plt.plot(x, y, ".")
plt.plot(x_grid, p)
plt.show()

检查模型参数：

print(plm.state_dict())
> OrderedDict([('breaks', tensor([[6.0033]])),
               ('linear.weight', tensor([[ 1.9999, 11.8892]])),
               ('linear.bias', tensor([2.9963]))])

神经网络的预测等价于：

def f(x): 
   return 1.9999*x + 11.8892*(x - 6.0033)*(x > 6.0033) + 2.9963

【参考方案11】：

您正在寻找线性树。它们是以广义和自动化的方式应用分段线性拟合（也适用于多变量和分类上下文）的最佳方法。

线性树与决策树不同，因为它们计算线性近似（而不是常数近似），在叶子中拟合简单的线性模型。

对于我的一个项目，我开发了linear-tree：一个 python 库，用于在叶子上构建带有线性模型的模型树。

linear-tree 被开发为可与 scikit-learn 完全集成。

from sklearn.linear_model import *
from lineartree import LinearTreeRegressor, LinearTreeClassifier

# REGRESSION
regr = LinearTreeRegressor(base_estimator=LinearRegression())
regr.fit(X, y)

# CLASSIFICATION
clf = LinearTreeClassifier(base_estimator=RidgeClassifier())
clf.fit(X, y)

LinearTreeRegressor 和 LinearTreeClassifier 以 scikit-learn BaseEstimator 的形式提供。它们是根据来自sklearn.linear_model 的线性估计器的数据构建决策树的包装器。 sklearn.linear_model 中可用的所有模型都可以用作线性估计器。

比较决策树和线性树：

考虑到您的数据，概括非常简单：

from sklearn.linear_model import LinearRegression
from lineartree import LinearTreeRegressor
import numpy as np
import matplotlib.pyplot as plt

X = np.array(
    [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15]
    ).reshape(-1,1)
y = np.array(
    [5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03]
    )

model = LinearTreeRegressor(base_estimator=LinearRegression())
model.fit(X, y)

plt.plot(X, y, ".", label='TRUE')
plt.plot(X, model.predict(X), label='PRED')
plt.legend()

【参考方案12】：

piecewise-regressionpython package 正好解决了这个问题。

import numpy as np
import matplotlib.pyplot as plt
import piecewise_regression

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15])
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03])

pw_fit = piecewise_regression.Fit(x, y, n_breakpoints=1)
pw_fit.plot()
plt.xlabel("x")
plt.ylabel("y")
plt.show()

它还给出了拟合的结果：

pw_fit.summary()

它通过实现 Muggeo 的迭代算法来工作。 More code examples here

带有一些噪音的示例。举个更有趣的例子，我们可以给 y 数据添加一些噪声并再次拟合：

y += np.random.normal(size=len(y)) * 5

pw_fit = piecewise_regression.Fit(x, y, n_breakpoints=1)
pw_fit.plot()