文章目录
多项式拟合非线性拟合多元拟合
Python科学计算：数组数据生成微积分插值
多项式拟合
所谓数据拟合，就是用一个系数待定的函数表达式，尽可能地逼近给定的一组数据。多项式拟合，顾名思义，即给定的函数是多项式。泰勒公式表明，任意连续函数都可以通过多项式来进行逼近，所以多项式的拟合能力并不局限于多项式函数。比如，正弦函数可以表示为
sin
x
=
x
−
x
3
3
!
+
x
5
5
!
−
⋯
\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\cdots
sinx=x−3!x3+5!x5−⋯
下面用一个5次多项式，来拟合一小段正弦函数，拟合结果如下图所示。
import numpy as np
x = np.arange(10)
y = np.sin(x)
abc = np.polyfit(x,y,5)
print(abc)
# [-7.293e-04, 3.183e-03, 1.169-01, -9.844e-01, 1.900, -4.222e-02]
fit = np.poly1d(abc)
X = np.linspace(0,9,100)
Y = fit(X)
plt.scatter(x, y, label='y=sin(x)')
plt.plot(X, Y, label='polyfit')
plt.legend()
plt.show()
在上述代码中，主要用到两个陌生的函数。
【polyfit】是numpy中用于多项式拟合的函数，其输入为
x
,
y
x,y
x,y以及最高次数deg，示例如下，由于在拟合时输入的最高项次数为5，故而其返回值是
x
5
x^5
x5的系数，根据Taylor公式计算可得
1
5
!
=
0.0083
\frac{1}{5!}=0.0083
5!1=0.0083，可见将Taylor公式截断到5次，与真实结果还是有一定出入的。
【poly1d】函数用于生成一个多项式函数，其输入的参数为多项式的系数。fit就是poly1d使用拟合结果生成的多项式函数，将
X
X
X代入其中，即可得到输出值。
非线性拟合
scipy作为专业的科学计算模块，如果不能突破numpy的限制，那就没什么意思了。其【optimize】模块中，提供了一个非线性拟合函数【curve_fit】，可以实现对自定义的函数表达式进行拟合。以高斯函数
y
=
a
exp
−
(
x
−
b
c
)
2
y = a\exp-(\frac{x-b}{c})^2
y=aexp−(cx−b)2w为例，其拟合结果
代码为
from scipy.optimize import curve_fit
def gauss(x, a, b, c):
return a*np.exp(-(x-b)**2/c**2)
x = np.arange(100)/10
y = gauss(x, 2, 5, 3) + np.random.rand(100)/10
abc, para = curve_fit(gauss, x, y)
print(abc)
# [2.03042233 5.01182397 3.10994351]
plt.scatter(x, y, marker='.', label="gauss")
plt.plot(x, gauss(x, *abc), lw=1, label="curve_fit")
plt.grid()
plt.legend()
plt.show()
【curve_fit】在调用时输入了三个参数，分别是拟合函数、自变量、因变量。返回值abc和para分别为拟合参数和拟合的协方差，最终得到拟合结果与预设值
2
,
0.5
,
3
2,0.5, 3
2,0.5,3十分接近。
除了这三个必选参数之外，【curve_fit】还提供了下列参数，可以进一步定制拟合流程。
p0 拟合参数初始值sigma 相对精度要求absolute_sigma 绝对精度要求check_finite 有限性检测开关bounds 拟合范围method 拟合方法，可选‘lm’, ‘trf’, ‘dogbox’，与least_squares函数中定义相同jac 雅可比矩阵，与least_squares中定义相同
多元拟合
尽管curve_fit的参数列表中，只给出了自变量xdata和因变量ydata这两组参数用于拟合，但并未规定xdata必须是一维数组，换言之，curve_fit具备多元函数拟合的潜力。
从其内部逻辑出发，对形如
y
=
f
(
x
)
y=f(x)
y=f(x)的函数在进行拟合时，必然会对比当前迭代的拟合值
y
^
\hat y
y^和输入数据
y
y
y的比对，即
δ
=
∣
y
^
−
y
∣
\delta=\vert\hat y-y\vert
δ=∣y^−y∣，而这步运算成立的关键，就是要求
y
^
\hat y
y^和
y
y
y的维度一致。在curve_fit中，要求二者必须是一维数组。
而
y
^
\hat y
y^是在迭代中得到的拟合函数函数
f
^
(
x
)
\hat f(x)
f^(x)计算得来，因此，如果想实现多元函数的拟合，至少要求
f
(
x
)
f(x)
f(x)的返回值是一维数组。
以二元高斯函数
z
=
a
exp
−
(
x
0
−
b
)
2
+
(
x
1
−
d
)
2
2
c
2
z=a\exp-\frac{(x_0-b)^2+(x_1-d)^2}{2c^2}
z=aexp−2c2(x0−b)2+(x1−d)2为例，在创建函数时，需要写成如下形式
def g2(x, a, b, c, d):
d2 = (x[0] - b) ** 2 + (x[1] - d) ** 2
r = a * np.exp(-d2 / (2 * c ** 2))
return r.ravel()
对此函数进行数值拟合，结果如下图所示。
代码如下和部分输出如下，可以发现拟合结果与预设的
a
b
c
d
abcd
abcd十分接近，也验证了这种方法的合理性。
# 生成原始数据
xx = np.indices([10, 10])
z = g2(xx, 10, 5, 2, 5) + np.random.normal(size=100)/100
# 数据拟合
abcd, para = curve_fit(g2, xx, z)
print(abcd)
# [10.00258587 5.00146314 1.99952885 5.00138184]
# 绘图
z = z.reshape(10, 10)
Z = g2(xx, *abcd).reshape(10,10)
ax = plt.subplot(projection='3d')
ax.scatter3D(xx[0], xx[1], z, color='red')
ax.plot_surface(xx[0], xx[1], Z, cmap='rainbow')
plt.show()
参考阅读
您阅读本篇文章共花了：
发表评论