最小二乘法曲线拟合原理及maab实现

1、曲线拟合(curve-fitting ):工程实践中，用测量到的一些离散的数据(X,yj,i求一个近似的函数(x)来拟合这组数据，要求所得的拟合曲线能最好的反映数据的基本趋势(即使(X)最好地逼近fx,而不必满足插值原则。因此没必 要取(x) =yi,只要使i (Xi) yi尽可能地小)。原理：给定数据点( Xi,yi)J0 ,1 , 2,. m o求近似曲线(x)。并且使得近似曲线与f x的偏差最小。(Xi ) yj ,i=1o常见的曲线拟合方法：1. 使偏差绝对值之和最小2. 使偏差绝对值最大的最小3. 使偏差平方和最小最小二乘法：按偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟

2、合曲线的方法，称为最小二乘法。 推导过程：1. 设拟合多项式为：2. 各点到这条曲线的距离之和，即偏差平方和如下：3问题转化为求待定系数ao.ak对等式右边求q偏导数，因而我们得到了：4、把这些等式化简并表示成矩阵的形式，就可以得到下面的矩阵：5. 将这个范德蒙得矩阵化简后可得到:6. 也就是说X*A=Y ,那么A = (X,*X)-1*X,*Y，便得到了系数矩阵A,同时，我们也就得到了拟合 曲线。MATLAB 实现：MATLAB提供了 polyfit ()函数命令进行最小二乘曲线拟合。调用格式：p=polyfit(x,y5n)p,s= polyfit(x,y,n)p,s,mu=polyfit

3、(x,y,n)x,y为数据点，n为多项式阶数，返回p为幕次从高到低的多项式系数向量po x 必须是单调的。矩阵s包括R (对x进行QR分解的三角元素)、df(自由度)、 normr(残差)用于生成预测值的误差估计。p,s,mu=polyfit(x,y,n)在拟合过程中，首先对x进行数据标准化处理，以在拟合中 消除量纲等影响，mu包含标准化处理过程中使用的x的均值和标准差。polyval() 为多项式曲线求值函数，调用格式:y=polyval(p,x)y,DELTA=polyval(p,x,s)y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。y,DELTA=polyval

4、(p,x,s)使用polyfit函数的选项输出s得出误差估计丫 DELTAo它 假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。贝ij Y DELTA将 至少包含50%的预测值。如下给定数据的拟合曲线：x=0.5,1.0,1.5,2.0,2.5,3.0,y=1.75,2.45,3.81,4.80,7.00,8.60o 解：MATLAB 程序如下： x=0.5,1.0,1.552.0,2.553.0;y=1.75,2.45,3.81,4.80,7.00,8.60;p=polyfit(x,y,2) x1 =0.5:0.05:3.0;y1=polyval(p,x1); plot(x

5、,y/*r,x15y1 ；-b*)运行结果如图 1 计算结果为：p =0.5614 0.8287 1.1560 即所得多项式为 y=0.5614xA2+0.08287x+1.15560图1最小二乘法曲线拟合示例对比检验拟合的有效 性：例：在0,n区间上对正弦函数进行拟合，然后在0,2 n区间画出图形，比较拟 合区间和非拟合区间的图形，考察拟合的有效性。在MATLAB中输入如下代码：clearx=0:0.1:pi; y=sin(x);p,mu=polyfit(x,y,9) x1=0:0.1:2*pi;y1=sin(x1);%实际曲线y2=polyval(p,x1);%根据由区间0到pi上进行拟合

6、得到的多项式计算0到2pi上的函 数值，%需要注意的是polyval ()返回的函数值在pi到2pi上并 没有 进行拟合plot(x1,y2,k*,x1,y1,'kJ)运行结果:P =0.0000 0.0000 -0.0003 0.0002 0.0080 0.0002 -0.1668 0.0000 1.0000 0.0000mu =R: 10x10 doubledf: 22 normr: 1.6178e-07MATLAB的最优化工具箱还提供了 Isqcurvefit ()函数命令进行最小二乘曲线 拟合 (Solve non li near curve-fitti ng (data-fi

7、tti ng) problems in least-squares sense)o 调用格式：x = lsqcurvefit(fu n5x0,xdata,ydata)x = lsqcurvefit(fu n,xO,xdata,ydata,lb,ub)x = lsqcurvefit(fu n,xO,xdata,ydata,lb,ub,optio ns)x = Isqcurvefit(problem)x3res norm = lsqcurvefit(.)x,res no rm5residual = lsqcurvefit(.)x,res norm,residual,exitflag = lsqcu

8、rvefit(.)x,res norm,residual,exitflag,output = Isqcurvefit()x,res norm,residual,exitflag,output,lambda=Isqcurvefit()xjes no rm,residual3exitflag,output,lambdajacobia n = xO 为初女台角军向量;xdata , yd ata 为;茜 足关系ydata=F(x, xdata)的数据；lb、ub为解向量的下界和上界，若没有指定界，贝Ilb=, ub=；options为指定的优化参数；fun 为拟合函数，其定义方式为:x = lsqc

9、urvefit(myfun,xO,xdata,ydata),其中 myfun 已定义为 function F =myfun(x,xdata)F =%计算x处拟合函数值fun的用法与前面相同；resnorm=sum (fun(x,xdata)-ydataF2),即在 x 处残差的平方和；residual=fun(x,xdata)-ydata ,即在 x 处 的残差；exitflag为终止迭代的条件； output为输出的优化信息；lambda为解x处的Lagrange乘子；jacobian为解x处拟合函数fun的jacobian矩阵。 例:lsqcurvefit()优化程序Data =.0.00

10、000.10005.89553.56390.20000.30002.51731.97900.40000.50000.60000.70000.80001.89901.39381.13591.00961.03430.90001.00000.84350.68561.10001.20001.30000.61000.53920.39461.40001.50001.60001.70000.39030.54740.34590.13701.80001.90002.00000.22110.17040.2636;t = Data(:,1);y = Data(:,2);% axis(0 2 -0.5 6)plot(

11、t,y；ro')title('Data points')% We would like to fit the function y = c(1)*exp(-lam(1)*t) + c(2)*exp(-la m(2)*t) to the data %The Isqcurvefit function solves this type of problem easily.%To begin, define the parameters in terms of one variable x:%x(1) = c(1)%x(2) = lam(1)%x(3) = c(2)%x(4)

12、= lam(2)%Then define the curve as a function of the parameters x and the data t:F = (x,xdata)x(1 )*exp(-x(2)*xdata) + x(3)*exp(-x(4)*xdata);xO = 1 1 1 0;x,res norm,exitflag,output = lsqcurvefit(F,xO,t5y)hold onplot(t,F(x,t)hold offFsumsquares = (x)sum(F(x,t) y).A2);opts = optimsetCLargeScale'/of

13、f);xunc3ressquared,eflag,outputu=.fminun c(Fsumsquares,xO,opts)fprintf(,There were %d iterations using fminunc；.'and %d using IsqcurvefitA outputu.iteratio ns,output, iteratio ns)fprintf(,There were %d function evaluations using fminunc,'and %d using Isqcurvefit.1, .outputu.f uncCoun t,outpu

14、t.f uncCount)type fitvectorx02 = 1 0;F2 = (x,t) fitvector(x,t,y);x2,resnorm2,exitflag2,output2 = lsqcurvefit(F2,xO2,t,y)fprintf(,There were %d function evaluations using the 2-d 1.'formulation, and %d using the 4-d formulation.1,. output2.funcCount5output.funcCount) xObad = 511 0;xbad,res no rmb

15、ad5exitflagbad5outputbad=.lsqcurvefit(F,xObad,t,y)hold onplot(t,F(xbad,t),'g,)legend('Data','Global fit','Bad local fit','Location*,'NE')hold offfpri ntf(,The residual norm at the good ending point is %f；.'and the residual norm at the bad ending point is %fJ,.res no rm,res no rmbad)displayE ndOfDemoMessage(mfile name)拟合效果如下：直线的最小二乘拟合:y= a+bx式中有两个待定参数，a代表截距，b代表斜率。对于等精度测量所得到的N组 数据(xi, yi), i= 1, 2,N, xi值被认为是准确的，所有的误差只联 系着yi。下 面利用最小二乘法把观测数据拟合为直线。用最小二乘法估计参数时，要求观测值yi的偏差的加权平方和为最小。对 于等精度观测值的直线拟合来说，可使下式的值最小:上式分别对a、b求偏导得：整理后得到方程组：解上述方程组便可求得直线参数a和