数学建模

 

编程知识点

模型库-scipy-拟合-自定义函数曲线拟合(leastsq)

 

Wang Haihua

🚅 🚋😜 🚑 🚔

1. 数据拟合

在科学研究和工程应用中经常通过测量、采样、实验等方法获得各种数据。对一组已知数据点集，通过调整拟合函数（曲线）的参数，使该函数与已知数据点集相吻合，这个过程称为数据拟合，又称曲线拟合。

插值和拟合都是根据一组已知数据点，求变化规律和特征相似的近似曲线的过程。但是插值要求近似曲线完全经过所有的给定数据点，而拟合只要求近似曲线在整体上尽可能接近数据点，并反映数据的变化规律和发展趋势。因此插值可以看作是一种特殊的拟合，是要求误差函数为 0 的拟合。

2. 自定义函数曲线拟合

2.1 scipy.optimize.curve_fit() 函数说明

• curve_fit() 使用非线性最小二乘法将自定义的拟合函数拟合到观测数据，不仅可以用于直线、二次曲线、三次曲线的拟合，而且可以适用于任意形式的自定义函数的拟合，使用非常方便。curve_fit() 允许进行单变量或多变量的自定义函数拟合。

• scipy.optimize.curve_fit(f,xdata,ydata,p0=None,sigma=None,absolute_sigma=False,check_finite=True,bounds=(-inf,inf),method=None,jac=None,kwargs)

主要参数：

• f：可调用的函数，自定义的拟合函数，具有一个或多个待定参数。拟合函数的形式为 func(x,p1,p2,...)，其中参数必须按照 (x,p1,p2,...) 的顺序排列，p1, p2,... 是标量不能表达为数组。
• xdata：n*m数组，n 为观测数据长度，m为变量个数。观测数据 xdata 可以是一维数组（单变量问题），也可以是多维数组（多变量问题）。
• ydata：数组，长度为观测数据长度 n。
• p0：可选项，待定参数 [p1,p2,...] 的初值，默认值无。

返回值：

• popt：待定参数的最小二乘估计值。
• pcov：参数的估计值 popt 的协方差，其对角线是各参数的方差。

2.2 ：单变量自定义函数曲线拟合

程序说明：

不同于 leastsq() 定义观测值与拟合函数值的误差函数，scipy.optimize.curve_fit() 直接定义一个自定义的拟合函数，更为直观和便于理解。 curve_fit() 定义一个拟合函数，函数名可以任意定义，但拟合函数的参数必须按照 (x,p1,p2,...) 的顺序排列，不能改变次序。p1, p2,... 是标量，不能写成数组。注意 leastsq() 中误差函数的参数必须按照 (p,x,y) 的顺序排列，与 curve_fit() 不同。 leastsq() 也可以对自定义的拟合函数进行最小二乘拟合。 由于本例程中自定义拟合函数使用了观测数据的实际模型，而不是通用的多项式函数或样条函数，因此拟合结果不仅能很好的拟合观测数据，而且能更准确地反映实际模型的趋势。 Python 例程：

# 6. 自定义函数曲线拟合：单变量
import numpy as np
import matplotlib.pyplot as plt  # 导入 Matplotlib 工具包
from scipy.optimize import leastsq, curve_fit  # 导入 scipy 中的曲线拟合工具

def fitfunc6(x, p0, p1, p2, p3):  # 定义拟合函数为自定义函数
    # p0, p1, p2, p3 = p  # 拟合函数的参数
    y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)
    return y

# def error6(p, x, y):  # 定义观测值与拟合函数值的误差函数
#     p0, p1, p2, p3 = p
#     err = fitfunc6(x, p0, p1, p2, p3) - y  # 计算残差
#     return err

# 创建给定数据点集 (x, yObs)
p0, p1, p2, p3 = [1.0, 1.2, 0.5, 0.8]  # y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)
x = np.linspace(-1, 2, 30)
y = fitfunc6(x, p0, p1, p2, p3)
np.random.seed(1)
yObs = y + 0.5*np.random.randn(x.shape[-1])  # 生成带有噪声的观测数据

# # 用 scipy.optimize.leastsq() 进行函数拟合
# pIni = [1, 1, 1, 1]  # 设置拟合函数的参数初值
# pFit, info = leastsq(error6, pIni, args=(x, yObs))  # 最小二乘法求拟合参数
# print("Data fitting of custom function by leastsq")
# print("y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)")
# print("p[0] = {:.4f}\np[1] = {:.4f}\np[2] = {:.4f}\np[3] = {:.4f}"
#       .format(pFit[0], pFit[1], pFit[2], pFit[3]))

# 用 scipy.optimize.curve_fit() 进行自定义函数拟合(单变量)
pFit, pcov = curve_fit(fitfunc6, x, yObs)
print("Data fitting of custom function by curve_fit:")
print("y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)")
print("p[0] = {:.4f}\np[1] = {:.4f}\np[2] = {:.4f}\np[3] = {:.4f}"
      .format(pFit[0], pFit[1], pFit[2], pFit[3]))
print("estimated covariancepcov:\n",pcov)
# 由拟合函数 fitfunc 计算拟合曲线在数据点的函数值
yFit = fitfunc6(x, pFit[0], pFit[1], pFit[2], pFit[3])

# 绘图
fig, ax = plt.subplots(figsize=(8,6))
ax.set_title("Data fitting of custom function")
plt.scatter(x, yObs, label="observed data")
plt.plot(x, y, 'r--', label="theoretical curve")
plt.plot(x, yFit, 'b-', label="fitting curve")
plt.legend(loc="best")
plt.show()

Data fitting of custom function by curve_fit:
y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)
p[0] = 0.9460
p[1] = 1.1465
p[2] = 0.8291
p[3] = 0.6008
estimated covariancepcov:
 [[ 0.01341654  0.00523061 -0.01645431  0.00455901]
 [ 0.00523061  0.02648836 -0.04442234  0.02821206]
 [-0.01645431 -0.04442234  0.20326672 -0.07482843]
 [ 0.00455901  0.02821206 -0.07482843  0.0388316 ]]

2.3 Python 例程：多变量自定义函数曲线拟合

程序说明：

scipy.optimize.curve_fit() 既可以用于单变量也可以用于多变量问题，本例程求解一个二元非线性拟合问题。 curve_fit() 定义一个拟合函数 fitfunc7(X, p0, p1, p2, p3)，函数名可以任意定义，但拟合函数的参数必须按照 (x,p1,p2,...) 的顺序排列，不能改变次序。p1, p2,... 是标量，不能写成数组。 curve_fit(fitfunc7, X, yObs) 中的 X 是 (n,m) 数组，n 是观测数据点集的长度，m 是变量个数。 Python 例程：

# 7. 自定义函数曲线拟合：多变量
import numpy as np
import matplotlib.pyplot as plt  # 导入 Matplotlib 工具包
from scipy.optimize import curve_fit  # 导入 scipy 中的曲线拟合工具

def fitfunc7(X, p0, p1, p2, p3):  # 定义多变量拟合函数, X 是向量
    # p0, p1, p2, p3 = p  # 拟合函数的参数
    y = p0 + p1*X[0,:] + p2*X[1,:] + p3*np.sin(X[0,:]+X[1,:]+X[0,:]**2+X[1,:]**2)
    return y

# 创建给定数据点集 (x,yObs)
p = [1.0, 0.5, -0.5, 5.0]  # 自定义函数的参数
p0, p1, p2, p3 = p  # y = p0 + p1*x1 + p2*x2 + p3*np.sin(x1+x2+x1^2+x2^2)
np.random.seed(1)
x1 = 2.0 * np.random.rand(8)  # 生成随机数组，长度为 8
x2 = 3.0 * np.random.rand(5)  # 生成随机数组，取值范围 (0,3.0)
xmesh1, xmesh2 = np.meshgrid(x1, x2)  # 生成网格点的坐标 xx,yy (二维数组)
xx1= xmesh1.reshape(xmesh1.shape[0]*xmesh1.shape[1], )  # 将网格点展平为一维数组
xx2= xmesh2.reshape(xmesh2.shape[0]*xmesh2.shape[1], )  # 将网格点展平为一维数组
X = np.vstack((xx1,xx2))  # 生成多变量数组，行数为变量个数
y = fitfunc7(X, p0, p1, p2, p3)  # 理论计算值 y=f(X,p)
yObs = y + 0.2*np.random.randn(y.shape[-1])  # 生成带有噪声的观测数据
print(x1.shape,x2.shape,xmesh1.shape,xx1.shape,X.shape)

# 用 scipy.optimize.curve_fit() 进行自定义函数拟合(多变量)
pFit, pcov = curve_fit(fitfunc7, X, yObs)  # 非线性最小二乘法曲线拟合
print("Data fitting of multivariable custom function")
print("y = p0 + p1*x1 + p2*x2 + p3*np.sin(x1+x2+x1^2+x2^2)")
for i in range(4):
    print("p[{:d}] = {:.4f}\tp[{:d}]_fit = {:.4f}".format(i, p[i], i, pFit[i]))

# 由拟合函数 fitfunc 计算拟合曲线在数据点的函数值
yFit = fitfunc7(X, pFit[0], pFit[1], pFit[2], pFit[3])

(8,) (5,) (5, 8) (40,) (2, 40)
Data fitting of multivariable custom function
y = p0 + p1*x1 + p2*x2 + p3*np.sin(x1+x2+x1^2+x2^2)
p[0] = 1.0000	p[0]_fit = 1.1316
p[1] = 0.5000	p[1]_fit = 0.5020
p[2] = -0.5000	p[2]_fit = -0.5906
p[3] = 5.0000	p[3]_fit = 5.0061

参考文献

• youcans Python小白的数学建模课 https://www.zhihu.com/column/c_1381900867424075776