Numpy 使用Scipy curve_fit进行分段函数拟合

在本文中，我们将介绍使用Numpy和Scipy库进行分段函数拟合的方法。分段函数拟合是一种数据拟合技术，使用分段函数可以更好地逼近非连续、不光滑或其他具有显著特征的数据集。我们将使用Scipy中的curve_fit函数来拟合具有斜率或曲率不连续变化的分段函数。

阅读更多：Numpy 教程

什么是分段函数？

分段函数是指在不同区间内使用不同的函数表达式的函数。这些函数表达式被称为分段函数的“分段”。例如，一个分段线性函数可以被写成如下形式：

$f(x) = \begin{cases} a_1 x + b_1&x \leq x_1\\a_2 x + b_2&x > x_1 \end{cases}$

其中， $a_1, b_1, a_2, b_2$ 是要拟合的系数， $x_1$ 是分段点。

分段函数可以更好地逼近一些非线性或不连续的数据。例如下图中的脉冲信号就可以用分段线性函数进行拟合。

使用curve_fit拟合分段函数

使用Python进行分段函数拟合的常用工具是Scipy库中的curve_fit函数。该函数需要传入一个自定义的拟合函数和需要拟合的数据。

下面的例子中，我们将使用以下分段二次函数来拟合一组带有噪声的数据：

$f(x) = \begin{cases} a_1 x^2 + b_1 x + c_1&x \leq x_1\\a_2 x^2 + b_2 x + c_2&x \leq x_2\\a_3 x^2 + b_3 x + c_3&x > x_2 \end{cases}$

由于该函数具有三个不同区间，我们需要传入三个不同区间的系数。

import numpy as np
from scipy.optimize import curve_fit

# 定义分段函数
def piecewise_quad(x, x0, x1, a0, b0, c0, a1, b1, c1, a2, b2, c2):
    # 多个分段函数的定义看起来有些冗长，但是只需要针对具体问题进行修改即可
    return np.piecewise(x, [x <= x0, (x > x0) & (x <= x1), x > x1],
                        [lambda x:a0*x**2 + b0*x + c0,
                        lambda x:a1*x**2 + b1*x + c1,
                        lambda x:a2*x**2 + b2*x + c2])

# 输入数据
xdata = np.linspace(0, 2.5, 200)
y = piecewise_quad(xdata, 0.8, 1.6, 2.5, -6, 5, 4, -6, 1, 4, 2, 3)

# 添加噪声并拟合
np.random.seed(2021)
y_noise = 0.2 * np.random.normal(size=xdata.size)
ydata = y + y_noise
popt, pcov = curve_fit(piecewise_quad, xdata, ydata)

# 绘制拟合结果
import matplotlib.pyplot as plt
plt.figure()
plt.plot(xdata, ydata, 'b-', label='data')
plt.plot(xdata, piecewise_quad(xdata, *popt), 'r-', label='fit')
plt.legend()
plt.show()