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

1、曲线拟合curve-fitting:工程实践中,用测量到的一些离散的数据x,yi,i0,1,2,.m求一个近似的函数x来拟合这组数据,要求所得的拟合曲线能最好的反映数据的根本趋势即使x最好地逼近fx,而不必满足插值原那么.因此没必要取xi=y只要使i国y尽可能地小.原理：给定数据点xi,yi0,1,2,.m.求近似曲线x.并且使得近似曲线与fx的偏差最小.近似曲线在该点处的偏差ixyi,i=1,2,.,m.常见的曲线拟合方法：1 .使偏差绝对值之和最小而“£|包|二2y三一匕口修2 .使偏差绝对值最大的最小minmax|=|那么»乂|中i3 .使偏差平方和最小mm'

2、r=i=l按偏差平方和最小的原那么选取拟合曲线,并且采取二项式方程为拟合曲线的方法,称为最小二乘法.推导过程:1 .设拟合多项式为:kxa0a1x.akx2 .各点到这条曲线的距离之和,即偏差平方和如下:及-三X民-M口+口1%+-+At.3 .问题转化为求待定系数ao.ak对等式右边求ai偏导数,因而我们得到了:昨-2%卜一.口x+.;班,01 =IJT-z£卜勤主如1+.上鼻/卜二o百库4、把这些等式化简并表示成矩阵的形式,就可以得到下面的矩阵:5 .将这个范德蒙得矩阵化简后可得到6 .也就是说X*A=Y,那么A=(X'*X)-1*X'*Y,便得到了系数矩阵A,同

3、时,我们也就得到了拟合曲线.MATLAB实现：MATLAB提供了polyfit()函数命令进行最小二乘曲线拟合.调用格式：p=polyfit(x,y,n)p,s=polyfit(x,y,n)p,s,mu=polyfit(x,y,n)x,y为数据点,n为多项式阶数,返回p为幕次从高到低的多项式系数向量p.x必须是单调的.矩阵s包括R6xx进彳TQR分解的三角元素)、df(自由度)、normr(残差)用于生成预测值的误差估计.p,s,mu=poly巾t(x,y,n)在拟合过程中,首先对x进行数据标准化处理,以在拟合中消除量纲等影响,mu包含标准化处理过程中使用的x的均值和标准差.polyval()

4、为多项式曲线求值函数,调用格式：y=polyval(p,x)y,DELTA=polyval(p,x,s)y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值.y,DELTA=polyval(p,x,s瘦用polyfit函数的选项输出s得出误差估计YDELTA它假设polyfit函数数据输入的误差是独立正态的,并且方差为常数.那么YDELTA!至少包含50%的预测值.如下给定数据的拟合曲线：x=,y=,.解：MATLAB程序如下:y=,;p=polyfit(x,y,2)x1=:;y1=polyval(p,x1);plot(x,y,'*r',x1,y1,'

5、-b')运行结果如图1计算结果为：p=即所得多项式为y=A2+9876543210.511.522.53图1最小二乘法曲线拟合例如比照检验拟合的有效性：例：在0,何区间上对正弦函数进行拟合,然后在0,2句区间画出图形,比拟拟合区间和非拟合区间的图形,考察拟合的有效性.在MATLAB中输入如下代码:clearx=0:pi;y=sin(x);p,mu=polyfit(x,y,9)x1=0:2*pi;y1=sin(x1);%际曲线y2=polyval(p,x1);%根据由区间0到pi上进行拟合得到的多项式计算0到2pi上的函数值,%需要注意的是polyval()返回的函数值在pi到2pi上并

6、没有进行拟合plot(x1,y2,'k*',x1,y1,'k-')运行结果：p=mu=R:10x10doubledf:22normr:0.80.60.40.2土-0.2-0.4-0.6-0.8-156*MATLAB的最优化工具箱还提供了Isqcurve巾t()函数命令进行最小二乘曲线拟合(Solvenonlinearcurve-fitting(data-fitting)problemsinleast-squaressense).调用格式：x=lsqcurvefit(fun,x0,xdata,ydata)x=lsqcurvefit(fun,x0,xdata,yda

7、ta,lb,ub)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)x=lsqcurvefit(problem)x,resnorm=lsqcurvefit(.)x,resnorm,residual=lsqcurvefit(.)x,resnorm,residual,exitflag=lsqcurvefit(.)x,resnorm,residual,exitflag,output=lsqcurvefit(.)x,resnorm,residual,exitflag,output,lambda=lsqcurvefit(.)x,resnorm,residual

8、,exitflag,output,lambda,jacobian=x0为初始解向量；xdata,ydata为满足关系ydata=F(x,xdata)的数据；lb、ub为解向量的下界和上界,假设没有指定界,那么lb=,ub=;options为指定的优化参数；fun为拟合函数,其定义方式为:x=lsqcurvefit(myfun,x0,xdata,ydata),其中myfun已定义为functionF=myfun(x,xdata)F=计算x处拟合函数值fun的用法与前面相同；resnorm=sum(fun(x,xdata)-ydata).A2),即在x处残差的平方和；residual=fun(x,

9、xdata)-ydata,即在x处的残差;exitflag为终止迭代的条件；output为输出的优化信息；lambda为解x处白LLagrange乘子；jacobian为解x处拟合函数fun的jacobian矩阵.例：lsqcurve巾t()优化程序Data=.;t=Data(:,1);y=Data(:,2);%axis(026)plot(t,y,'ro')title('Datapoints')%Wewouldliketofitthefunctiony=c(1)*exp(-lam(1)*t)+c(2)*exp(-lam(2)*t)tothedata%Thelsq

10、curvefitfunctionsolvesthistypeofproblemeasily.%Tobegin,definetheparametersintermsofonevariablex:%x(1)=c(1)%x(2)=lam(1)%x(3)=c(2)%x(4)=lam(2)%Thendefinethecurveasafunctionoftheparametersxandthedatat:F=(x,xdata)x(1)*exp(-x(2)*xdata)+x(3)*exp(-x(4)*xdata);x0=1110;x,resnorm,-,exitflag,output=lsqcurvefit

11、(F,x0,t,y)holdonplot(t,F(x,t)holdoffFsumsquares=(x)sum(F(x,t)-y).A2);opts=optimset('LargeScale','off);xunc,ressquared,eflag,outputu=.fminunc(Fsumsquares,x0,opts)fprintf('Therewere%diterationsusingfminunc,'.'and%dusinglsqcurvefit.n',.,fprintf('Therewere%dfunctionevalua

12、tionsusingfminunc,'.'and%dusinglsqcurvefit.',.,typefitvectorx02=10;F2=(x,t)fitvector(x,t,y);x2,resnorm2,exitflag2,output2=lsqcurvefit(F2,x02,t,y)fprintf('Therewere%dfunctionevaluationsusingthe2-d'.'formulation,and%dusingthe4-dformulation.',.,x0bad=5110;xbad,resnormbad,-,e

13、xitflagbad,outputbad=.lsqcurvefit(F,x0bad,t,y)holdonplot(t,F(xbad,t),'g')legend('Data','Globalfit','Badlocalfit','Location','NE')holdofffprintf('Theresidualnormatthegoodendingpointis%f,''andtheresidualnormatthebadendingpointis%f.',.resn

14、orm,resnormbad)displayEndOfDemoMessage(mfilename)拟合效果如下:Datapoints54321DataGlobalfitBadlocalfit-00.20.40.60.811.21.41.61.82直线的最小二乘拟合:y=a+bx式中有两个待定参数,a代表截距,b代表斜率.对于等精度测量所得到的N组数据xi,yi,i=1,2,N,xi值被认为是准确的,所有的误差只联系着yi.下面利用最小二乘法把观测数据拟合为直线.用最小二乘法估计参数时,要求观测值yi的偏差的加权平方和为最小.对于等精度观测值的直线拟合来说,可使下式的值最小：£乂-.+%1=1上式分别对a、b求偏导得：6-Vr2y£2乂-.+如=2之凹-口-如=0,诙IMX却N而Z氏一5+如=一22稣一Q+如卜产01整理后得到方程组：小毛二&#