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

曲线拟合（curve-fitting）：工程实践中，用测量到的一些离散的数据},...2,1,0),,{(miyxii求一个近似的函数)(x来拟合这组数据，要求所得的拟合曲线能最好的反映数据的基本趋势（即使)(x最好地逼近xf，而不必满足插值原则。因此没必要取)(ix=iy，只要使iiiyx)(尽可能地小）。原理：给定数据点},...2,1,0),,{(miyxii。求近似曲线)(x。并且使得近似曲线与xf的偏差最小。近似曲线在该点处的偏差iiiyx)(，i=1,2,...,m。常见的曲线拟合方法:1.使偏差绝对值之和最小2.使偏差绝对值最大的最小3.使偏差平方和最小最小二乘法：按偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟合曲线的方法,称为最小二乘法。推导过程：1.设拟合多项式为：kkxaxaax...)(102.各点到这条曲线的距离之和，即偏差平方和如下：3.问题转化为求待定系数0a...ka对等式右边求ia偏导数，因而我们得到了：.......4、把这些等式化简并表示成矩阵的形式，就可以得到下面的矩阵：5.将这个范德蒙得矩阵化简后可得到:6.也就是说X*A=Y，那么A=(X'*X)-1*X'*Y，便得到了系数矩阵A，同时，我们也就得到了拟合曲线。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包括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(p,x,s)使用polyfit函数的选项输出s得出误差估计YDELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。则YDELTA将至少包含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.60]。解：MATLAB程序如下：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.60];p=polyfit(x,y,2)x1=0.5:0.05:3.0;y1=polyval(p,x1);plot(x,y,'*r',x1,y1,'-b')运行结果如图1计算结果为：p=0.56140.82871.1560即所得多项式为y=0.5614x^2+0.08287x+1.155600.511.522.53123456789图1最小二乘法曲线拟合示例对比检验拟合的有效性：例：在[0,π]区间上对正弦函数进行拟合，然后在[0,2π]区间画出图形，比较拟合区间和非拟合区间的图形，考察拟合的有效性。在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上进行拟合得到的多项式计算0到2pi上的函数值，%需要注意的是polyval（）返回的函数值在pi到2pi上并没有进行拟合plot(x1,y2,'k*',x1,y1,'k-')运行结果：p=0.00000.0000-0.00030.00020.00800.0002-0.16680.00001.00000.0000mu=R:[10x10double]df:22normr:1.6178e-0701234567-1-0.8-0.6-0.4-0.200.20.40.60.81MATLAB的最优化工具箱还提供了lsqcurvefit（）函数命令进行最小二乘曲线拟合(Solvenonlinearcurve-fitting(data-fitting)problemsinleast-squaressense)。调用格式：x=lsqcurvefit(fun,x0,xdata,ydata)x=lsqcurvefit(fun,x0,xdata,ydata,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,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).^2)，即在x处残差的平方和;residual=fun(x,xdata)-ydata，即在x处的残差;exitflag为终止迭代的条件;output为输出的优化信息;lambda为解x处的Lagrange乘子;jacobian为解x处拟合函数fun的jacobian矩阵。例：lsqcurvefit()优化程序Data=...[0.00005.89550.10003.56390.20002.51730.30001.97900.40001.89900.50001.39380.60001.13590.70001.00960.80001.03430.90000.84351.00000.68561.10000.61001.20000.53921.30000.39461.40000.39031.50000.54741.60000.34591.70000.13701.80000.22111.90000.17042.00000.2636];t=Data(:,1);y=Data(:,2);%axis([02-0.56])plot(t,y,'ro')title('Datapoints')%Wewouldliketofitthefunctiony=c(1)*exp(-lam(1)*t)+c(2)*exp(-lam(2)*t)tothedata%Thelsqcurvefitfunctionsolvesthistypeofproblemeasily.%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(F,x0,t,y)holdonplot(t,F(x,t))holdoffFsumsquares=@(x)sum((F(x,t)-y).^2);opts=optimset('LargeScale','off');[xunc,ressquared,eflag,outputu]=...fminunc(Fsumsquares,x0,opts)fprintf(['Therewere%diterationsusingfminunc,'...'and%dusinglsqcurvefit.\n'],...outputu.iterations,output.iterations)fprintf(['Therewere%dfunctionevaluationsusingfminunc,'...'and%dusinglsqcurvefit.'],...outputu.funcCount,output.funcCount)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.'],...output2.funcCount,output.funcCount)x0bad=[5110];[xbad,resnormbad,~,exitflagbad,outputbad]=...lsqcurvefit(F,x0bad,t,y)holdonplot(t,F(xbad,t),'g')legend('Data','Globalfit','Badlocalfit','Location','NE')holdofffprintf(['Theresidualnormatthegoodendingpointis%f,'...'andtheresidualnormatthebadendingpointis%f.'],...resnorm,resnormbad)displayEndOfDemoMessage(mfilename)拟合效果如下：00.20.40.60.811.21.41.61.820123456DatapointsDataGlobalfitBadlocalfit直线的最小二乘拟合：y＝a+bx式中有两个待定参数，a代表截距，b代表斜率。对于等精度测量所得到的N组数据（xi，yi），i＝1，2……，N，xi值被认为是准确的，所有的误差只联系着yi。下面利用最小二乘法把观测数据拟合为直线。用最小二乘法估计参数时，要求观测值yi的偏差的加权平方和为最小。对于等精度观测值的直线拟合来说，可使下式的值最小：上式分别对a、b求偏导得：整理后得到方程组：解上述方程组便可求得直线参数a和b的最佳估计值。1、可看成是一阶多项式拟合，跟前面曲线拟合方法一样。2、利用linefit()函数进行最小二乘的直线拟合使用：clearx=[0.511.522.53];y=[1.752.453.814.888.6];[k,b]=linefit(x,y)%得到斜率k和常数by1=polyval([k,b],x);plot(x,y1,’k-’,x,y,’k*’)MATLAB一元到多元线性回归方程的计算和检验_百度文库=_71b8kDsSUfnYSXWKsvL_X2DAWqojz3