数学建模
编程知识点
模型库-scipy-拟合-非线性函数拟合
Wang Haihua
🚅 🚋😜 🚑 🚔
1. 数据拟合¶
在科学研究和工程应用中经常通过测量、采样、实验等方法获得各种数据。对一组已知数据点集，通过调整拟合函数（曲线）的参数，使该函数与已知数据点集相吻合，这个过程称为数据拟合，又称曲线拟合。
插值和拟合都是根据一组已知数据点，求变化规律和特征相似的近似曲线的过程。但是插值要求近似曲线完全经过所有的给定数据点，而拟合只要求近似曲线在整体上尽可能接近数据点，并反映数据的变化规律和发展趋势。因此插值可以看作是一种特殊的拟合，是要求误差函数为 0 的拟合。
以下介绍 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)
主要参数：
a：(m,n) 数组，表示方程 Ax=b 的左侧。b：(m,) 数组，表示方程 Ax=b 的右侧。
返回值：
x：最小二乘的解，指多项式函数的系数。 注意：lstsq() 函数中求解方程 Ax=b，A 是指由观测数据 按多项式阶数转换为矩阵,而 x 是指多项式函数的系数，详见例程。
2.2 Python 例程：指数函数拟合¶
程序说明：
scipy.optimize.leastsq()本质上是求解带有待定参数的误差函数最小化问题，因此可以用于指数函数的最小二乘拟合。类似地，原理上 leastsq() 也可以用于其它形式非线性函数的拟合问题，但对于某些函数拟合误差可能会比较大。leastsq()以子函数error3(p, x, y)来定义观测值与拟合函数值的误差函数，以动态参数 args 的方式传递观测数据 (x, yObs) 。leastsq()中误差函数的函数名可以任意定义，但误差函数的参数必须按照(p,x1,x2,y)的顺序排列，不能改变次序。
In [8]:
# 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)
In [9]:
# 由给定数据点集 (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)
Data fitting of exponential function y = p0 + p1 * np.exp(-p2*x) p[0] = 0.5216 p[1] = 2.5742 p[2] = 1.6875
In [11]:
# 绘图
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")
Out[11]:
<matplotlib.legend.Legend at 0x1becf5810a0>
In [12]:
plt.show()
2.3 案例3（多项式函数拟合）¶
程序说明：
scipy.optimize.leastsq()本质上是求解带有待定参数的误差函数最小化问题，因此可以用于多项式函数的最小二乘拟合，使用方法与线性拟合、指数拟合类似。由于
leastsq()要以子函数 error(p, x, y) 来定义观测值与拟合函数值的误差函数，在比较不同阶数的多项式函数拟合时需要定义多个对应的误差函数，比较繁琐。scipy.linalg.lstsq()只要传入观测数据 (x,yObs)，并将 x 按多项式阶数转换为 X，即可求出多项式函数的系数，不需要定义拟合函数或误差函数，非常适合比较不同阶数的多项式函数拟合的效果。对于相同阶数的多项式函数，leastsq() 与 lstsq() 的参数估计和拟合结果是相同的。
- 增大多项式的阶数，可以减小拟合曲线与观测数据的误差平方和，但也更容易导致过拟合，虽然能更好地拟合训练数据，但并不能真实反映数据的总体规律，因而对于训练数据以外的测试数据的拟合效果反而降低了。
In [15]:
# 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
In [16]:
# 创建给定数据点集 (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]))
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
In [ ]:
# 用 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))
In [17]:
plt.legend(loc="best")
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
Out[17]:
<matplotlib.legend.Legend at 0x1becf6485b0>
In [18]:
plt.show()
2.4 案例4（样条曲线拟合）¶
程序说明：
scipy.interpolate.UnivariateSpline()类是一种基于固定数据点创建函数的方法，使用样条曲线拟合到给定的数据点集。UnivariateSpline类由已知数据点集生成样条插值函数 y=spl(x)，通过调用样条插值函数可以计算指定 x 的函数值 f(x)。- `UnivariateSpline 类既可以进行数据插值，也可以进行拟合。参数 s=0 表示数据插值，样条曲线必须通过所有数据点；s>0 表示数据拟合，默认 s= len(w)。 通过 set_smoothing_factor(sf) 设置光滑因子，可以对样条拟合函数进行调节，使拟合曲线更好地反映观测数据特征，避免过拟合。 Python 例程：
In [23]:
# 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]) # 生成带有噪声的观测数据
In [24]:
# 由给定数据点集 (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")
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]
Out[24]:
<matplotlib.legend.Legend at 0x1becf6c51c0>
In [25]:
plt.show()
参考文献
- youcans Python小白的数学建模课 https://www.zhihu.com/column/c_1381900867424075776