Python科学计算：数组💯数据生成💯微积分💯插值

多项式拟合

所谓数据拟合，就是用一个系数待定的函数表达式，尽可能地逼近给定的一组数据。多项式拟合，顾名思义，即给定的函数是多项式。泰勒公式表明，任意连续函数都可以通过多项式来进行逼近，所以多项式的拟合能力并不局限于多项式函数。比如，正弦函数可以表示为

$$\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\cdots$$

下面用一个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$以及最高次数deg，示例如下，由于在拟合时输入的最高项次数为5，故而其返回值是$x^5$的系数，根据Taylor公式计算可得$\frac{1}{5!}=0.0083$，可见将Taylor公式截断到5次，与真实结果还是有一定出入的。

【poly1d】函数用于生成一个多项式函数，其输入的参数为多项式的系数。fit就是poly1d使用拟合结果生成的多项式函数，将$X$代入其中，即可得到输出值。

非线性拟合

scipy作为专业的科学计算模块，如果不能突破numpy的限制，那就没什么意思了。其【optimize】模块中，提供了一个非线性拟合函数【curve_fit】，可以实现对自定义的函数表达式进行拟合。以高斯函数$y=a\exp -(\frac{x-b}{c}{)}^2$w为例，其拟合结果

代码为

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$十分接近。

除了这三个必选参数之外，【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)$的函数在进行拟合时，必然会对比当前迭代的拟合值$\hat{y}$和输入数据$y$的比对，即$\delta =|\hat{y}-y|$，而这步运算成立的关键，就是要求$\hat{y}$和$y$的维度一致。在curve_fit中，要求二者必须是一维数组。

而$\hat{y}$是在迭代中得到的拟合函数函数$\hat{f}(x)$计算得来，因此，如果想实现多元函数的拟合，至少要求$f(x)$的返回值是一维数组。

以二元高斯函数$z=a\exp -\frac{(x_0-b{)}^2+(x_1-d{)}^2}{2c^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()

对此函数进行数值拟合，结果如下图所示。

代码如下和部分输出如下，可以发现拟合结果与预设的$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()