最小二乘法曲线拟合原理及m a t l a b实现Modified by JEEP on December 26th, 2020.曲线拟合(curve-fitting ):工程实践中,用测量到的一些离散的数据},...2,1,0),,{(m i y x i i =求一个近似的函数)(x ϕ来拟合这组数据,要求所得的拟合曲线能最好的反映数据的基本趋势(即使)(x ϕ最好地逼近()x f ,而不必满足插值原则。
因此没必要取)(i x ϕ=i y ,只要使i i i y x -=)(ϕδ尽可能地小)。
原理:给定数据点},...2,1,0),,{(m i y x i i =。
求近似曲线)(x ϕ。
并且使得近似曲线与()x f 的偏差最小。
近似曲线在该点处的偏差i i i y x -=)(ϕδ,i=1,2,...,m 。
常见的曲线拟合方法:1.使偏差绝对值之和最小2.使偏差绝对值最大的最小3.使偏差平方和最小最小二乘法:按偏差平方和最小的原则选取拟合曲线,并且采取二项式方程为拟合曲线的方法,称为最小二乘法。
推导过程:1. 设拟合多项式为:2. 各点到这条曲线的距离之和,即偏差平方和如下:3. 问题转化为求待定系数0a ...k a 对等式右边求i a 偏导数,因而我们得到了: .......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得出误差估计Y DELTA。
它假设polyfit函数数据输入的误差是独立正态的,并且方差为常数。
则Y DELTA将至少包含50%的预测值。
如下给定数据的拟合曲线:x=[,,,,,],y=[,,,,,]。
解:MATLAB程序如下:x=[,,,,,];y=[,,,,,];p=polyfit(x,y,2)x1=::;y1=polyval(p,x1);plot(x,y,'*r',x1,y1,'-b')运行结果如图1计算结果为:p =即所得多项式为y=^2++图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上并没有进行拟合plot(x1,y2,'k*',x1,y1,'k-')运行结果:p =mu =R: [10x10 double]df: 22normr:MATLAB的最优化工具箱还提供了lsqcurvefit()函数命令进行最小二乘曲线拟合(Solve nonlinear curve-fitting (data-fitting) problems in least-squares sense)。
调用格式: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已定义为 function F = 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 = ...[];t = Data(:,1);y = Data(:,2);% axis([0 2 6])plot(t,y,'ro')title('Data points')%We would like to fit the function y = c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t) to the data %The lsqcurvefit 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) = 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);x0 = [1 1 1 0];[x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,t,y)hold onplot(t,F(x,t))hold offFsumsquares = @(x)sum((F(x,t) - y).^2);opts = optimset('LargeScale','off');[xunc,ressquared,eflag,outputu] = ...fminunc(Fsumsquares,x0,opts)fprintf(['There were %d iterations using fminunc,' ...' and %d using lsqcurvefit.\n'], ...,fprintf(['There were %d function evaluations using fminunc,' ...' and %d using lsqcurvefit.'], ...,type fitvectorx02 = [1 0];F2 = @(x,t) fitvector(x,t,y);[x2,resnorm2,~,exitflag2,output2] = lsqcurvefit(F2,x02,t,y)fprintf(['There were %d function evaluations using the 2-d ' ...'formulation, and %d using the 4-d formulation.'], ...,x0bad = [5 1 1 0];[xbad,resnormbad,~,exitflagbad,outputbad] = ...lsqcurvefit(F,x0bad,t,y)hold onplot(t,F(xbad,t),'g')legend('Data','Global fit','Bad local fit','Location','NE')hold offfprintf(['The residual norm at the good ending point is %f,' ...' and the residual norm at the bad ending point is %f.'], ...resnorm,resnormbad)displayEndOfDemoMessage(mfilename)拟合效果如下:直线的最小二乘拟合:y=a+bx式中有两个待定参数,a代表截距,b代表斜率。
对于等精度测量所得到的N组数据(xi,yi),i=1,2……,N,xi值被认为是准确的,所有的误差只联系着yi。