Numpy分段线性拟合与n个断点

在本文中，我们将介绍numpy库中分段线性拟合的相关知识以及如何使用它进行n个断点的分段线性拟合。

阅读更多：Numpy 教程

什么是分段线性拟合？

分段线性拟合是一种在多个区间内使用不同的线性方程来拟合数据的方法。在这个方法中，先选定n个断点，再在这些断点上进行线性拟合，从而生成n+1段线性曲线。

假设我们有如下数据：

x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
y = [1, 2, 2.5, 3, 4, 5, 7, 9, 10, 12]

我们希望找到一条最优的线性曲线来拟合这些数据。但是如果直接使用线性曲线，很难拟合曲线中的突变点。因此我们可以使用分段线性拟合来解决这个问题。

使用numpy进行分段线性拟合

numpy提供了一个方便的函数可以用于分段线性拟合，即numpy.piecewise()函数。该函数可以接受一个数组、一个或多个条件和每个条件对应的函数，从而生成分段线性曲线。

比如，以下代码可以生成一个有两个断点的线性曲线：

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(10)
y = np.array([1, 2, 2.5, 3, 4, 5, 7, 9, 10, 12])

# 定义两个断点
b1, b2 = 3, 7

# 定义对应的函数
f1 = lambda x: x * 0 + np.mean(y[:b1])
f2 = lambda x: (x - b1) * (y[b2] - y[b1]) / (b2 - b1) + y[b1]
f3 = lambda x: x * 0 + np.mean(y[b2:])

# 使用piecewise进行分段线性拟合
y_fit = np.piecewise(x, [x<b1, (x>=b1)&(x<b2), x>=b2], [f1, f2, f3])

# 绘制原始数据和拟合曲线
plt.plot(x, y, 'o')
plt.plot(x, y_fit, '-')
plt.show()

可以看到，分段线性曲线在突变点处比线性曲线更准确地拟合了数据。

如何确定断点的位置？

在分段线性拟合中，确定断点的位置是非常重要的。通常情况下，我们可以尝试多个断点的位置，从而寻找最合适的断点。

手动确定断点

我们可以通过直观观察，手动确定很多个断点位置，从而寻找一个较为合适的拟合函数。

使用手动选点的方法，可以根据对数据的了解和对图像的感性认识，更好地确定断点位置。

基于信息标准确定断点

我们可以使用贝叶斯信息标准（BIC）确定最优的断点数量。BIC是基于最大似然估计的一个信息准则，可以在不过拟合数据的同时减少模型中的参数数量。

对于我们想要拟合的数据，设可能的断点数量为n，设n个断点位置为b1, b2, …, bn。然后，我们可以使用numpy.piecewise()函数进行拟合，并计算出拟合误差的BIC值。

以下是一个示例代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress

x = np.arange(10)
y = np.array([1, 2, 2.5, 3, 4, 5, 7, 9, 10, 12])

def piecewise(x, b):
    """使用n个断点拟合数据"""
    y_fit = np.zeros(len(x))
    for i in range(len(y_fit)):
        if x[i] < b[0]:
            y_fit[i] = np.mean(y[:b[0]])
        elif x[i] >= b[-1]:
            y_fit[i] = np.mean(y[b[-1]:])
        else:
            j = np.where(x[i] < b)[0][0]
            y_fit[i] = (x[i] - b[j-1]) * (y[b[j]] - y[b[j-1]]) / (b[j] - b[j-1]) + y[b[j-1]]
    return y_fit

def bic(y, y_fit, k):
    """计算BIC值"""
    n = len(y)
    mse = np.sum((y - y_fit)**2) / n
    bic = n * np.log(mse) + k * np.log(n)
    return bic

# 选择断点的数量
n = 5

# 生成所有可能的断点位置
b_list = []
for i in range(2, len(x)-2):
    for j in range(i+1, len(x)-1):
        for k in range(j+1, len(x)):
            b_list.append([i, j, k])

# 定义一个最小BIC的值和对应的断点位置
bic_min = np.inf
b_min = []

# 遍历所有的断点位置，选择最优的拟合函数
for b in b_list:
    y_fit = piecewise(x, b)
    k = 2 * len(b)  # 计算参数数量
    bic_val = bic(y, y_fit, k)
    if bic_val < bic_min:
        bic_min = bic_val
        b_min = b

# 绘制原始数据和拟合曲线
y_fit = piecewise(x, b_min)
plt.plot(x, y, 'o')
plt.plot(x, y_fit, '-')
plt.show()

print("最佳拟合的断点位置为：", b_min)

这段代码会输出最优的断点位置，以及对应的拟合曲线。根据计算结果，最优的断点位置为[3, 5, 7]

由于分段线性拟合是一种非常灵活的方法，因此我们可以根据实际需要对断点位置进行调整，从而得到更好的拟合效果。