在科学研究和工程应用中经常通过测量、采样、实验等方法获得各种数据。对一组已知数据点集，通过调整拟合函数（曲线）的参数，使该函数与已知数据点集相吻合，这个过程称为数据拟合，又称曲线拟合。
插值和拟合都是根据一组已知数据点，求变化规律和特征相似的近似曲线的过程。但是插值要求近似曲线完全经过所有的给定数据点，而拟合只要求近似曲线在整体上尽可能接近数据点，并反映数据的变化规律和发展趋势。因此插值可以看作是一种特殊的拟合，是要求误差函数为 0 的拟合。
数据拟合问题，可以从不同角度进行分类：
数据拟合通过调整拟合函数中的待定参数，从整体上接近已知的数据点集。
这是一个优化问题，决策变量是拟合函数的待定参数，优化目标是观测数据与拟合函数的函数值之间的某种误差指标。典型的优化目标是拟合函数值与观测值的误差平方和；当观测数据的重要性不同或分布不均匀时，也可以使用加权误差平方和作为优化目标。
数据拟合的基本方法是最小二乘法。对于观测数据 ( x i , y i ) , i = 1 , . . n，将观测值 y i与拟合函数 y = f ( x , p ) 的计算值 f ( x i ) 的误差平方和最小作为优化问题的目标函数：
( p 1 , ⋯ p m ) 是拟合函数中的待定参数。
对于线性拟合问题，设拟合函数为直线 f ( x ) = p 0 + p 1 ∗ x ， 由极值的必要条件 $ partial J/partial p_j = 0,; (j=0,1)$ 可以解出系数 p 0 , p 1 ：
对于多变量线性最小二乘问题，设拟合函数为直线 f ( x ) = p 0 + p 1 ∗ x 1 + ⋯ + p m ∗ x m ， 类似地，可以解出系数 p 0 , p 1 , ⋯ p m 。
对于非线性函数的拟合问题，通常也是按照最小二乘法的思路，求解上述误差平方和最小化这个非线性优化问题，常用的具体算法有搜索算法和迭代算法两类。
数据拟合是常用算法，Python 语言的很多工具包都提供了数据拟合方法，常用的如 Scipy、Numpy、Statsmodel、Scikit-learn 工具包都带有数据拟合的函数与应用。
Scipy 是最常用的 Python 工具包，本系列中非线性规划、插值方法也都是使用 Scipy 工具包实现，因此仍以 Scipy 工具包讲解数据拟合问题。
Scipy 工具包对于不同类型的数据拟合问题，提供了不同的函数或类。由于 Scipy 工具包是多个团队合作完成，而且经过了不断更新，因此调用不同函数和方法时的函数定义和参数设置有所差异，往往使小白感到困惑。
本文对单变量、多变量线性最小二乘拟合，指数函数、多项式函数、样条函数的非线性拟合，单变量、多变量的自定义函数拟合问题进行分析、给出完整例程和结果，数据拟合从此无忧。
线性最小二乘拟合是最简单和最常用的拟合方法。scipy.optimize 工具箱中的 leastsq()、lsq_linear()，scipy.stats 工具箱中的 linregress()，都可以实现线性最小二乘拟合。
leastsq() 根据观测数据进行最小二乘拟合计算，只需要观测值与拟合函数值的误差函数和待定参数 的初值，返回拟合函数中的待定参数 ( p 1 , ⋯ p m ) ，但不能提供参数估计的统计信息。leastsq() 可以进行单变量或多变量线性最小二乘拟合，对变量进行预处理后也可以进行多项式函数拟合。
scipy.optimize.leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)
主要参数：
返回值：
linregress() 根据两组观测数据 (x,y) 进行线性最小二乘回归，不仅返回拟合函数中的待定参数 ( p 1 , p 1 )，而且可以提供参数估计的各种统计信息，但只能进行单变量线性拟合。
scipy.stats.linregress(x, y=None, alternative=‘two-sided’)
主要参数：
返回值：
程序说明：
Python 例程：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 1. 单变量线性拟合：最小二乘法 scipy.optimize.leastsq
import numpy as np
import matplotlib.pyplot as plt # 导入 Matplotlib 工具包
from scipy.optimize import leastsq # 导入 scipy 中的最小二乘法拟合工具
from scipy.stats import linregress # 导入 scipy 中的线性回归工具
def fitfunc1(p, x): # 定义拟合函数为直线
p0, p1 = p # 拟合函数的参数
y = p0 + p1*x # 拟合函数的表达式
return y
def error1(p, x, y): # 定义观测值与拟合函数值的误差函数
err = fitfunc1(p,x) - y # 误差
return err
# 创建给定数据点集 (x,yObs)
p = [2.5, 1.5] # y = p[0] + p[1] * x
x = np.array([0., 0.5, 1.5, 2.5, 4.5, 5.5, 7.5, 8.0, 8.5, 9.0, 10.0])
y = p[0] + p[1] * x # 理论值 y
np.random.seed(1)
yObs = y + np.random.randn(x.shape[-1]) # 生成带有噪声的观测数据
# print(x.shape, y.shape, yObs.shape)
# 由给定数据点集 (x,y) 求拟合函数的参数 pFit
p0 = [1, 1] # 设置拟合函数的参数初值
pFit, info = leastsq(error1, p0, args=(x,yObs)) # 最小二乘法求拟合参数
print("Data fitting with Scipy.optimize.leastsq")
print("y = p[0] + p[1] * x")
print("p[0] = {:.4f}np[1] = {:.4f}".format(pFit[0], pFit[1]))
# 由拟合函数 fitfunc 计算拟合曲线在数据点的函数值
yFit = fitfunc1(pFit,x)
# 比较：线性回归，可以返回斜率，截距，r 值，p 值，标准误差
slope, intercept, r_value, p_value, std = linregress(x, yObs)
print("nLinear regress with Scipy.stats.linregress")
print("y = p[0] + p[1] * x")
print("p[0] = {:.4f}".format(intercept)) # 输出截距 intercept
print("p[1] = {:.4f}".format(slope)) # 输出斜率 slope
print("r^2_value: {:.4f}".format(r_value**2)) # 输出 r^2 值
print("p_value: {:.4f}".format(p_value)) # 输出 p 值
print("std: {:.4f}".format(std)) # 输出标准差 std
# 绘图
fig, ax = plt.subplots(figsize=(8,6))
ax.text(8,3,"youcans-xupt",color='gainsboro')
ax.set_title("Data fitting with linear least squares")
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 with Scipy.optimize.leastsq
y = p[0] + p[1] * x
p[0] = 2.2688
p[1] = 1.5528
Linear regress with Scipy.stats.linregress
y = p[0] + p[1] * x
p[0] = 2.2688
p[1] = 1.5528
r^2_value: 0.9521
p_value: 0.0000
std: 0.1161
程序说明：
Python 例程：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 2. 多变量线性拟合：最小二乘法 scipy.optimize.leastsq
import numpy as np
import matplotlib.pyplot as plt # 导入 Matplotlib 工具包
from scipy.optimize import leastsq # 导入 scipy 中的最小二乘法工具
def fitfunc2(p, x1, x2): # 定义拟合函数为直线
p0, p1, p2 = p # 拟合函数的参数
y = p0 + p1*x1 + p2*x2 # 拟合函数的表达式
return y
def error2(p, x1, x2, y): # 定义观测值与拟合函数值的误差函数
err = fitfunc2(p, x1, x2) - y # 计算残差
return err
# 创建给定数据点集 (x,yObs)
np.random.seed(1)
p = [2.5, 1.5, -0.5] # y = p[0] + p[1] * x1 + p[2] * x2
x1 = np.array([0., 0.5, 1.5, 2.5, 4.5, 5.5, 7.5, 8.0, 8.5, 9.0, 10.0])
x2 = np.array([0., 1.0, 1.5, 2.2, 4.8, 5.0, 6.3, 6.8, 7.1, 7.5, 8.0])
z = p[0] + p[1]*x1 + p[2]*x2 # 理论值 z
zObs = z + np.random.randn(x1.shape[-1]) # 生成带有噪声的观测数据
print(x1.shape, z.shape, zObs.shape)
# 由给定数据点集 (x,z) 求拟合函数的参数 pFit
p0 = [1, 1, 1] # 设置拟合函数的参数初值
pFit, info = leastsq(error2, p0, args=(x1,x2,zObs)) # 最小二乘法求拟合参数
print("Data fitting with Scipy.optimize.leastsq:")
print("z = p[0] + p[1]*x1 + p[1]*x2")
print("p[0]={:.4f}np[1]={:.4f}np[2]={:.4f}".format(pFit[0], pFit[1], pFit[2]))
# 由拟合函数 fitfunc 计算拟合曲线在数据点的函数值
zFit = fitfunc2(pFit, x1, x2)
# 多元线性回归：最小二乘法(OLS)
import statsmodels.api as sm
x0 = np.ones(x1.shape[-1]) # 截距列 x0=[1,...1]
X = np.column_stack((x0, x1, x2)) # (nSample,3): [x0, x1, x2]
model = sm.OLS(zObs, X) # 建立 OLS 模型: y = b0*x0 + b1*x1 + b2*x2 + e
results = model.fit() # 返回模型拟合结果
zFit = results.fittedvalues # 模型拟合的 y值
print(results.summary()) # 输出回归分析的摘要
print("nOLS model: y = b0*x0 + b1*x1 + b2*x2")
print('Parameters: ', results.params) # 输出：拟合模型的系数程序运行结果：
Data fitting with Scipy.optimize.leastsq
z = p[0] + p[1]*x1 + p[1]*x2
p[0]=2.6463
p[1]=2.2410
p[2]=-1.3710
OLS Regression Results
OLS model: y = b0*x0 + b1*x1 + b2*x2
Parameters: [ 2.64628055 2.24100973 -1.37104475]
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 2.6463 0.942 2.808 0.023 0.473 4.820
x1 2.2410 1.043 2.148 0.064 -0.165 4.647
x2 -1.3710 1.312 -1.045 0.326 -4.396 1.654
==============================================================================非线性函数是非常广泛的概念。本节讨论指数函数、多项式函数和样条函数三种常用的通用形式的非线性函数拟合问题，分别使用了 Scipy 工具包中的 scipy.optimize.leastsq()、scipy.linalg.lstsq() 和
scipy.interpolate.UnivariateSpline() 函数。
scipy.optimize.leastsq() 的使用方法已在本文 2.1 中进行了介绍，
scipy.interpolate.UnivariateSpline() 的使用方法在《22. 插值方法》文中进行了介绍，以下介绍 scipy.linalg.lstsq() 函数。
lstsq() 函数只要传入观测数据 (x,yObs)，并将 x 按多项式阶数转换为 X，即可求出多项式函数的系数，不需要定义拟合函数或误差函数，非常适合比较不同阶数的多项式函数拟合的效果。
scipy.linalg.lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False, check_finite=True, lapack_driver=None)
主要参数：
返回值：
注意：lstsq() 函数中求解方程 Ax=b，A 是指由观测数据 x i x_ix 按多项式阶数转换为矩阵 ( x i 0 , x i 1 , . . . x i m ) , i = 1 , ,n，b 是指 y i , i = 1 ,,n，而 x 是指多项式函数的系数，详见例程。
程序说明：
Python 例程：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 3. 非线性函数拟合：指数函数拟合(exponential function)
import numpy as np
import matplotlib.pyplot as plt # 导入 Matplotlib 工具包
from scipy.optimize import leastsq # 导入 scipy 中的最小二乘法工具
def fitfunc3(p, x): # 定义拟合函数
p0, p1, p2 = p # 拟合函数的参数
y = p0 + p1 * np.exp(-p2*x) # 拟合函数的表达式
return y
def error3(p, x, y): # 定义残差函数
err = fitfunc3(p,x) - y # 残差
return err
# 创建给定数据点集 (x,yObs)
p = [0.5, 2.5, 1.5] # y = p0 + p1 * np.exp(-p2*x)
x = np.linspace(0, 5, 50)
y = fitfunc3(p, x)
np.random.seed(1)
yObs = y + 0.2*np.random.randn(x.shape[-1]) # 生成带有噪声的观测数据
# print(x.shape, y.shape, yObs.shape)
# 由给定数据点集 (x,y) 求拟合函数的参数 pFit
p0 = [1, 1, 1] # 设置拟合函数的参数初值
pFit, info = leastsq(error3, p0, args=(x,yObs)) # 最小二乘法求拟合参数
print("Data fitting of exponential function")
print("y = p0 + p1 * np.exp(-p2*x)")
print("p[0] = {:.4f}np[1] = {:.4f}np[2] = {:.4f}".format(pFit[0], pFit[1], pFit[2]))
# 由拟合函数 fitfunc 计算拟合曲线在数据点的函数值
yFit = fitfunc3(pFit, x)
# 绘图
fig, ax = plt.subplots(figsize=(8,6))
ax.set_title("Data fitting of exponential 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 exponential function
y = p0 + p1 * np.exp(-p2*x)
p[0] = 0.5216
p[1] = 2.5742
p[2] = 1.6875
程序说明：
Python 例程：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 4. 非线性函数拟合：多项式函数拟合(Polynomial function)
import numpy as np
import matplotlib.pyplot as plt # 导入 Matplotlib 工具包
from scipy.linalg import lstsq # 导入 scipy 中的 linalg, stats 函数库
from scipy.optimize import leastsq # 导入 scipy 中的最小二乘法工具
def fitfunc4(p, x): # 定义拟合函数
p0, p1, p2, p3 = p # 拟合函数的参数
y = p0 + p1*x + p2*x*x + p3*x*x*x # 拟合函数的表达式
return y
def error4(p, x, y): # 定义观测值与拟合函数值的误差函数
err = fitfunc4(p,x) - y # 残差
return err
# 创建给定数据点集 (x,yObs)
p = [1.0, 1.2, 0.5, 0.8] # y = p0 + ((x*x-p1)**2+p2) * np.sin(x*p3)
func = lambda x: p[0]+((x*x-p[1])**2+p[2])*np.sin(x*p[3]) # 定义 y=f(x)
x = np.linspace(-1, 2, 30)
y = func(x) # 计算已知数据点的理论值 y =f(x)
np.random.seed(1)
yObs = y + 0.5*np.random.randn(x.shape[-1]) # 生成带有噪声的观测数据 yObs
# 绘图
fig, ax = plt.subplots(figsize=(8,6))
ax.set_title("Polynomial fitting with least squares")
plt.scatter(x, yObs, label="observed data")
plt.plot(x, y, 'c--', label="theoretical curve")
# 用 scipy.optimize.leastsq() 进行多项式函数拟合
p0 = [1, 1, 1, 1] # 设置拟合函数的参数初值
pFit, info = leastsq(error4, p0, args=(x,yObs)) # 最小二乘法求拟合参数
yFit = fitfunc4(pFit, x) # 由拟合函数 fitfunc 计算拟合曲线在数据点的函数值
ax.plot(x, yFit, '-', label='leastsq')
print("Polynomial fitting by scipy.optimize.leastsq")
print("y = p[0] + p[1]*x + p[2]*x^2 +p[3]*x^3") # 拟合函数的表达式
print("p[0] = {:.4f}np[1] = {:.4f}np[2] = {:.4f}np[3] = {:.4f}".format(pFit[0], pFit[1], pFit[2], pFit[3]))
# 用 scipy.linalg.lstsq() 进行多项式函数拟合
print("nPolynomial fitting by scipy.linalg.lstsq")
print("y = w[0] + w[1]*x + w[2]*x^2 + ... + w[m]*x^m") # 拟合函数的表达式
# 最小二乘法多项式数据拟合，求解多项式函数的系数 W=[w[0],...w[m]]
for order in range(1,5):
# order = 3 # 多项式阶数, m
X = np.array([[(xi ** i) for i in range(order + 1)] for xi in x])
Y = np.array(yObs).reshape((-1, 1))
W, res, rnk, s = lstsq(X, Y) # 最小二乘法求解 A*W = b
print("order={:d}".format(order))
for i in range(order+1): # W = [w[0],...w[order]]
print("tw[{:d}] = {:.4f}".format(i, W[i,0]))
# 由拟合函数 fitfunc 计算拟合曲线在数据点的函数值
yFit = X.dot(W) # y = w[0] + w[1]*x + w[2]*x^2 + ... + w[m]*x^m
# 绘图：n 次多项式函数拟合曲线
ax.plot(x, yFit, '-', label='order={}'.format(order))
plt.legend(loc="best")
plt.show()程序运行结果：
Polynomial fitting by scipy.optimize.leastsq
y = p[0] + p[1]*x + p[2]*x^2 +p[3]*x^3
p[0] = 0.9586
p[1] = -0.4745
p[2] = -0.3581
p[3] = 1.1721
Polynomial fitting by scipy.linalg.lstsq
y = w[0] + w[1]*x + w[2]*x^2 + ... + w[m]*x^m
order=1
w[0] = 1.0331
w[1] = 1.7354
order=2
w[0] = 0.2607
w[1] = 0.3354
w[2] = 1.4000
order=3
w[0] = 0.9586
w[1] = -0.4745
w[2] = -0.3581
w[3] = 1.1721
order=4
w[0] = 1.0131
w[1] = 1.6178
w[2] = -1.1039
w[3] = -1.5209
w[4] = 1.3465
程序说明：
scipy.interpolate.UnivariateSpline() 类是一种基于固定数据点创建函数的方法，使用样条曲线拟合到给定的数据点集。
UnivariateSpline 类由已知数据点集生成样条插值函数 y=spl(x)，通过调用样条插值函数可以计算指定 x 的函数值 f(x)。
UnivariateSpline 类既可以进行数据插值，也可以进行拟合。参数 s=0 表示数据插值，样条曲线必须通过所有数据点；s>0 表示数据拟合，默认 s= len(w)。
通过 set_smoothing_factor(sf) 设置光滑因子，可以对样条拟合函数进行调节，使拟合曲线更好地反映观测数据特征，避免过拟合。
Python 例程：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 5 非线性函数拟合：样条函数拟合(Spline function)
import numpy as np
import matplotlib.pyplot as plt # 导入 Matplotlib 工具包
from scipy.interpolate import UnivariateSpline # 导入 scipy 中的样条插值工具
# 创建给定数据点集 (x,yObs)
# y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3)
np.random.seed(1)
p0, p1, p2, p3 = [1.0, 1.2, 0.5, 0.8] # 拟合函数的参数
x = np.linspace(-1, 2, 30) # 生成已知数据点集的 x
y = p0 + ((x*x-p1)**2+ p2) * np.sin(x*p3) # 生成理论值 y
yObs = y + 0.5*np.random.randn(x.shape[-1]) # 生成带有噪声的观测数据
# 由给定数据点集 (x,y) 求拟合函数的参数 fSpl
fSpl = UnivariateSpline(x, yObs) # 三次样条插值，默认 s= len(w)
coeffs = fSpl.get_coeffs() # Return spline coefficients
print("Data fitting with spline function")
print("coeffs of 3rd spline function:n ", coeffs)
# 由拟合函数 fitfunc 计算拟合曲线在数据点的函数值
yFit = fSpl(x) # 由插值函数 fSpl1 计算插值点的函数值 yFit
# 对拟合函数 fitfunc 进行平滑处理
fSpl.set_smoothing_factor(10) # 设置光滑因子 sf
yS = fSpl(x) # 由插值函数 fSpl(sf=10) 计算插值点的函数值 ys
coeffs = fSpl.get_coeffs() # 平滑处理后的参数
print("coeffs of 3rd spline function (sf=10):n ", coeffs)
# 绘图
fig, ax = plt.subplots(figsize=(8,6))
ax.set_title("Data fitting with spline function")
plt.scatter(x, yObs, label="observed data")
plt.plot(x, y, 'r--', label="theoretical curve")
plt.plot(x, yFit, 'b-', label="3rd spline fitting")
plt.plot(x, yS, 'm-', label="smoothing spline")
plt.legend(loc="best")
plt.show()程序运行结果：
Data fitting with spline function
coeffs of 3rd spline function:
[-0.09707885 3.66083026 -4.20416235 7.95385344]
coeffs of 3rd spline function (sf=10):
[0.41218039 0.52795588 1.6248287 0.76540737 8.49462738]
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) **
主要参数：
返回值：
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 ]]结果分析：
程序说明：
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 例程：
# mathmodel25_v1.py
# Demo25 of mathematical modeling algorithm
# Demo of curve fitting with Scipy
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-03
# 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])程序运行结果：
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
estimated covariancepcov:
[[ 9.51937904e-03 -2.82863223e-03 -5.26393413e-03 -8.51457970e-04]
[-2.82863223e-03 4.88275894e-03 9.39281331e-05 3.73832161e-04]
[-5.26393413e-03 9.39281331e-05 3.86701646e-03 4.65766686e-04]
[-8.51457970e-04 3.73832161e-04 4.65766686e-04 1.85374067e-03]]
关注微信