前言拟合在数学建模中应用广泛,而最小二乘拟合是一种常用的拟合方法。数据拟合,又称曲线拟合,俗称拉曲线,是通过数学方法将已有数据代入数学表达式的一种方法。科学和工程问题可以通过采样、实验等方法获取一些离散数据。根据这些数据,我们往往希望得到与已知数据匹配的连续函数(即曲线)或更稠密的离散方程。这个过程叫做拟合。                                                                 --摘自百度百科最小二乘拟合,我们从原理出发,给出三个示例并结合python代码,从而加深对最小二乘拟合的理解。1.PrincipleFirstofall,leastsquaresfittingisanidea,通过最小二乘拟合、线性拟合、多项式拟合、指数拟合等可以进行任意形式的拟合。拟合:拟合用于解决为一系列已知的\(x_i,y_i\)求\(y=f(x)\)方程的问题。最小二分思路:我们通过拟合找到\(\phi(x)\),使得误差平方和Q=\(\sum(\phi(x_i)-y_i)^2\)最小,则为认为\(\phi(x)\)约等于我们需要的\(f(x)\)。最小化误差平方和的原理Q=\(\sum(\phi(x_i)-y_i)^2\)就是最小二分思想。\(\phi(x)\)的表达式:如何求\(\phi(x)\),可以选择\(\phi(x)=f(e^x),\phi(x)=f(x^n)+f(x^{n-1})+....\)等方程进行拟合。这时候会生成\(a*e^x+b\)中的a、b等一些参数。误差平方和Q:那么如何求Q的最小值呢?先介绍一下高等数学的内容$$\begin{cases}Q_a=0\\Q_b=0\&+函数连续\右箭头求极值点+边界点\右箭头极大点\\Q_c=0\\......\end{cases}$$极值的充分条件是Q对a,b,c...的偏导数等于0且Q类型连续。通过比较极值和边界值来确定最大值。这里求Q的最小值,通过求Q对a,b,c...等于0的偏导数,可以得到a,b,c...的参数,从而求Q的值。有n个方程和n个未知数,可以求解。因此,根据a,b,c...等参数的取值,可以得到\(\phi(x)\)的方程,用来画图比较拟合预测值和真正的价值。2.python代码——多项式拟合主要是np.polyfit和np.poly1d的使用。importnumpyasnpiimportmatplotlib.pyplotaspltfrompylabimportmmplmppl.rcParams["font.sans-serif"]=["SimHei"]#这两行解决字体问题x=np.arange(1,101)#更改datay=np.array([11,5,4,7,6,6,5,7,#改变数据13,6,5,7,12,5,4,6,9,5,5,11,29,21,17,20,27,13,9,10,16,6,5,7,11,5,5,6,12,7,7,10,15,10,9,11,15,10,10,16,26,21,23,36,50,45,45,49,57,43,40,44,52,43,42,45,52,41,39,41,48,35,34,35,42,34,36,43,55,48,54,65,80,70,74,85,101,89,88,90,100,87,88,89,104,89,89,90,106,96,94,99])z1=np.polyfit(x,y,deg=100)#deg=100,100次多项式,返回值为系数p1=np.poly1d(z1)#通过多项式系数,返回方程print(p1)#输出方程print(np.polyval(p1,12))#预测print(np.polyval(z1,13))#两种方法都可以,p1或者z1都可以作为参数,12、13为x,输入x得到预测值y_pred=p1(x)#预测值plt.plot(x,y,'*',label='原始值e')plt.plot(x,y_pred,'r',color='green',label='预测值')plt.title('多项式拟合')plt.xlabel('xlable')plt.ylabel('ylabel')plt.legend(loc=3,borderaxespad=0.,bbox_to_anchor=(0,0))#绘制图例plt.show()有兴趣的可以调整多项式的次数来研究过拟合的情况,这里让我解释更多。以下是多项式拟合情况——任意形式的拟合指数拟合。以前的拟合方法不能拟合其他形式的方程如\(f(e^x)\)。下面主要是通过plt.cur_fit来拟合。查找函数并记下函数返回类型。代码中的fun()函数可以设置任何我们想要的功能,从而实现任意形式的拟合。importnumpyasnpimportmathimportmatplotlib.pyplotaspltfromscipy.optimizeimportcurve_fitfrompylabimportmplmppl.rcParams["font.sans-serif"]=["SimHei"]#这两行解决字体问题deffun(x,a,b,c,d,e):#自定义任意类型函数returna*np.exp(x)+b*x+c#returna*b**x+c+d/np.exp(x)+e/x**5;if__name__=="__main__":x=np.arange(1,105)y=np.array([11,5,4,7,6,6,5,7,13,6,5,7,12,5,4,6,9,5,5,11,29,21,17,20,27,13,9,10,16,6,5,7,11,5,5,6,12,7,7,10,15,10,9,11,15,10,10,16,26,21,23,36,50,45,45,49,57,43,40,44,52,43,42,45,52,41,39,41,48,35,34,35,42,34,36,43,55,48,54,65,80,70,74,85,101,89,88,90,100,87,88,89,104,89,89,90,106,96,94,99,109,99,96,102])popt,pcov=curve_fit(fun,x,y)#popt返回值为残差最小时的参数,即最佳参数y_pred=[fun(i,popt[0],popt[1],popt[2],popt[3],popt[4])foriinx]#带入x和参数得到y的预测值print(popt)#输出参数plt.plot(x,y,'*',label='原始值')plt.plot(x,y_pred,'-',label='预测值',color='绿色')plt.xlabel('xlabel')plt.ylabel('ylabel')plt.title('指数拟合')plt.legend(loc=3,borderaxespad=0.,bbox_to_anchor=(0,0))plt.show()下面是指数函数拟合的情况,因为数据的问题差距有点大,但是明白了原理就不用关心这些细节了。替换下面的fun函数即可完成幂函数拟合。这里我顺便改了下数据。importnumpyasnpimportmathimportmatplotlib.pyplotaspltfromscipy.optimizeimportcurve_fitfrompylabimportmplmppl.rcParams["font.sans-serif"]=["SimHei"]#这两行解决字体问题deffun(x,a,b,c,d,e):#自定义任意类型函数returna/x+b#returna*np.exp(x)+b*x+c#returna*b**x+c+d/np.exp(x)+e/x**5;if__name__=="__main__":x=np.arange(1,6)y=np.array([1,0.78,0.61,0.55,0.42])popt,pcov=curve_fit(fun,x,y)#popt返回值为残差最小时的参数,即最佳参数y_pred=[fun(i,popt[0],popt[1],popt[2],popt[3],popt[4])foriinx]#带入x和参数得到y的预测值print(popt)#输出参数plt.plot(x,y,'*',label='原始值')plt.plot(x,y_pred,'-',label='预测值',color='green')plt.xlabel('xlabel')plt.ylabel('ylabel')plt.title('幂函数拟合')plt.legend(loc=3,borderaxespad=0.,bbox_to_anchor=(0,0))plt.show()如图:好了,到此结束~记录下我的先t写Blog,以后再接再厉~
