拟合
已知一组二维数据，即平面上$n$个点$(x_i,y_i)(i=1,2,…,n), x_i$互不相同，求一个函数(曲线)$y=f(x)$，使$f(x)$在某种准则下与所有数据最为接近，即为拟合，记
\(\delta_i=f(x_i)-y_i, i=1,2,...,n\) 称$ \delta_i$为拟合函数$f(x)$在$x_i $点处的偏差(参差).为了让数据最为接近，可以使用平方和最小作为判定规则
\(J=\sum^n_{i=1}(f(x_i)-y_i)^2\)
使J达到最小，这种方法称为最小二乘法
线性最小二乘拟合
非线性最小二乘拟合
拟合函数的选择
一般可以先画出散点图，然后判断用什么拟合函数
- 若数据分布接近直线，应使用线性函数$f(x)=a_1x+a_2$拟合
- 若数据分布接近抛物线，应使用二次多项式$f(x)=a_1x^2+a_2x+a_3$拟合
- 若数据开始上升快，然后变慢，应使用双曲线型函数$f(x)=\frac{x}{a_1x+a_2}$或指数型函数$f(x)=a_1e^{\frac{-a_2}{x}}$拟合
- 若数据开始下降快，然后变慢，可使用$f(x)=\frac{1}{a_1x+a_2}, f(x)=\frac{1}{a_1x^2+a_2} 或 f(x)=a_1e^{-a_2x}$等函数拟合
- 常使用的非线性函数有对数函数$y=a_1+a_2lnx$，S型曲线函数$y=\frac{1}{a+be^{-x}}$
实例
Numpy的ployfit函数可以做多项式函数的拟合，scipy.optimize模块的leastsq函数和curve_fit都可以拟合函数
ployfit
下表为待拟合的数据
| x_i | 0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 |
| :— | :— | :— | :— | :— | :— | :— | :— | :— | :— | :— | :— |
| y_i | -0.447 | 1.978 | 3.28 | 6.16 | 7.08 | 7.34 | 7.66 | 9.56 | 9.48 | 9.30 | 11.2 |
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
xdata = np.linspace(0,1,11)
ydata = np.array([-0.447, 1.978, 3.28, 6.16, 7.08, 7.34, 7.66, 9.56, 9.48, 9.30, 11.2])
p = np.polyfit(xdata, ydata, 2) # 拟合二次多项式，返回p为高次到低次系数
yhat = np.polyval(p, [0.25, 0.35]) # 使用p对列表中数据做预测
plt.plot(xdata, ydata, '*', xdata, np.polyval(p, xdata), '-')
plt.show()
curve_fit
curve_fit函数和leastsq函数差别并不大，顺便说下scipy有很多方法都集成在一个函数中，有些函数和这种集成的函数并没有什么区别(代码都一样，或者直接调用)
popt, pcov = curve_fit(func, xdata, ydata)
参数
- func: 要拟合的函数
- xdata, ydata: xy轴数据
返回值
- popt: 拟合的参数
- pocv: 参数的协方差矩阵
# 拟合二次多项式
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
%matplotlib inline
y = lambda x, a, b, c : a*x**2 + b*x + c
xdata = np.arange(0, 1.1, 0.1)
ydata = np.array([-0.447, 1.978, 3.28, 6.16, 7.08, 7.34, 7.66, 9.56, 9.48, 9.30, 11.2])
popt, pcov = curve_fit(y, xdata, ydata)
print('拟合的参数:', popt)
print('预测值:', y(np.array([0.25, 0.35]), *popt))
# ynew = lambda x: *popt[0]*x**2 + popt[1]*x + *popt[2]
plt.plot(xdata, ydata, '*', xdata, y(xdata, *popt))
plt.show()
可见两者差别不大。。。
- 使用表中数据拟合$z=ae^{bx}+cy^2$
| x | 6 | 2 | 6 | 7 | 4 | 2 | 5 | 9 |
| :— | :— | :— | :— | :— | :— | :— | :— | :— |
| y | 4 | 9 | 5 | 3 | 8 | 5 | 8 | 2 |
| z | 5 | 2 | 1 | 9 | 7 | 4 | 3 | 3 |
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from mpl_toolkits import mplot3d
from matplotlib.ticker import MaxNLocator
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
xdata = np.array([6, 2, 6, 7, 4, 2, 5, 9])
ydata = np.array([4, 9, 5, 3, 8, 5, 8, 2])
zdata = np.array([5, 2, 1, 9, 7, 4, 3, 3])
xydata = np.vstack((xdata, ydata))
def Pfun(t, a, b, c):
return a*np.exp(b*t[0]) + c*t[1]**2
popt, pocv = curve_fit(Pfun ,xydata, zdata)
print('a, b, c拟合值为:', popt)
X, Y = np.meshgrid(xdata, ydata)
Z = Pfun([X, Y], popt[0], popt[1], popt[2])
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet, linewidth=0)
fig.colorbar(surf)
title = ax.set_title("fitting results")
title.set_y(1.01)
ax.view_init(40, 50)
fig.tight_layout()
- 使用模拟数据拟合曲面$x=x=e^{\frac{-(x-\mu_1)^2+(y-\mu_2)^2}{2\delta^2}}$，其中$\mu_1=1, \mu_2=2, \delta=3$，并画出图形 ```python from mpl_toolkits import mplot3d import numpy as np from scipy.optimize import curve_fit import matplotlib.pyplot as plt
m = 200 n = 300 x = np.linspace(-6, 6, m) y = np.linspace(-8, 8, n) x2, y2 = np.meshgrid(x, y) x3 = np.reshape(x2, (1, -1)) y3 = np.reshape(y2, (1, -1)) xy = np.vstack((x3, y3))
def Pfun(t, m1, m2, s): return np.exp(-((t[0] - m1) ** 2 + (t[1] - m2) ** 2) / (2 * s ** 2)) z = Pfun(xy, 1, 2, 3) zr = z + 0.2 * np.random.normal(size=z.shape) # 为数据添加噪声
popt, pcov = curve_fit(Pfun, xy, zr) # 拟合参数 zn = Pfun(xy, *popt) zn2 = np.reshape(zn, x2.shape)
ax = plt.axes(projection='3d') # 三维对象 ax.plot_surface(x2, y2, zn2, cmap='gist_rainbow') ax.view_init(20, 50) plt.savefig(‘3dfitting') plt.show() ``` 