在数学建模过程中,通常要处理由试验、测量得到的大量数据或一些过于复杂而不便于计算的函数表达式,针对此情况,很自然的想法就是,构造一个简单的函数作为要考察数据或复杂函数的近似。插值和拟合就可以解决这样的问题。给定一组数据,需要确定满足特定要求的曲线(或曲面),如果所求曲线通过所给定有限个数据点,这就是插值。
得到简单实用的近似函数,这就是曲线拟合。插值和拟合都是根据一组数据构造一个函数作为近似。
插值
Lagrange插值
这式子显然就是拉格朗日式,名称由此而来。
#程序文件名Pfun7_1.py
def h(x,y,a):
s=0.0
for i in range(len(y)):
t=y[i]
for j in range(len(y)):
if i !=j:
t*=(a-x[j])/(x[i]-x[j])
s +=t
return s
用Python求解插值问题【第40页】
scipy.interpolate模块有一维插值函数interp1d,二维插值函数interp2d,多维插值函数interpn、interpnd。
interp1d的基本调用格式为
interp1d(x, y, kind='linear')
其中kind的取值是字符串,指明插值方法,kind的取值可以为:linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic'等,这里的'zero', 'slinear', 'quadratic' and 'cubic'分别指的是0阶、1阶、2阶和3阶样条插值。
拟合【第52页】
最小二乘的意思,偏差最小。判定公式:
线性最小二乘拟合 第54页
非线性最小二乘拟合 58页
拟合函数的选择
数据拟合时,首要也是最关键的一步就是选取恰当的拟合函数。如果能够根据问题的背景通过机理分析得到变量之间的函数关系,那么只需估计相应的参数即可。但很多情况下,问题的机理并不清楚。此时,一个较为自然的方法是先做出数据的散点图,从直观上判断应选用什么样的拟合函数。
一般来讲,如果数据分布接近于直线,则宜选用线性函数拟合;如果数据分布接近于抛物线,则宜选用二次多项式拟合;如果数据分布特点是开始上升较快随后逐渐变缓,则宜选用双曲线型函数或指数型函数,即用:
如果数据分布特点是开始下降较快随后逐渐变缓,则宜选用:
数据拟合的Python实现【62】
NumPy库中的多项式拟合函数polyfit;scipy.optimize模块中的函数leastsq,curve_fit都可以拟合函数。下面介绍polyfit和curve_fit的使用方法。
#程序文件名Pex7_7.py
from numpy import polyfit, polyval, array, arange
from matplotlib.pyplot import plot,show,rc
x0=arange(0, 1.1, 0.1)
y0=array([-0.447, 1.978, 3.28, 6.16, 7.08, 7.34, 7.66, 9.56, 9.48, 9.30, 11.2])
p=polyfit(x0, y0, 2) #拟合二次多项式
print("拟合二次多项式的从高次幂到低次幂系数分别为:",p)
yhat=polyval(p,[0.25, 0.35]); print("预测值分别为:", yhat)
rc('font',size=16)
plot(x0, y0, '*', x0, polyval(p, x0), '-'); show()
curve_fit的调用格式为:
popt, pcov = curve_fit(func, xdata, ydata)
其中func是拟合的函数,xdata是自变量的观测值,ydata是函数的观测值,返回值popt是拟合的参数,pcov是参数的协方差矩阵。
#程序文件名Pex7_8.py
import numpy as np
from scipy.optimize import curve_fit
y=lambda x, a, b, c: a*x**2+b*x+c
x0=np.arange(0, 1.1, 0.1)
y0=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, x0, y0)
print("拟合的参数值为:", popt)
print("预测值分别为:", y(np.array([0.25, 0.35]), *popt))
运行结果如下:
拟合的参数值为: [-9.81083901 20.12929291 -0.03167108]
预测值分别为: [4.38747471 5.81175366]
#程序文件Pex7_10.py
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) #拟合参数
print("三个参数的拟合值分别为:",popt)
zn=Pfun(xy, *popt) #计算拟合函数的值
zn2=np.reshape(zn, x2.shape)
plt.rc('font',size=16)
ax=plt.axes(projection='3d') #创建一个三维坐标轴对象
ax.plot_surface(x2, y2, zn2,cmap='gist_rainbow')
plt.savefig("figure7_10.png", dpi=500); plt.show()