最小二乘法曲线拟合原理及matlab实现.docx
《最小二乘法曲线拟合原理及matlab实现.docx》由会员分享,可在线阅读,更多相关《最小二乘法曲线拟合原理及matlab实现.docx(13页珍藏版)》请在冰豆网上搜索。
![最小二乘法曲线拟合原理及matlab实现.docx](https://file1.bdocx.com/fileroot1/2023-7/26/68d59ca4-1ae0-4ee8-aee3-19e187360665/68d59ca4-1ae0-4ee8-aee3-19e1873606651.gif)
最小二乘法曲线拟合原理及matlab实现
曲线拟合(curve-fitting):
工程实践中,用测量到的一些离散的数据
求一个近似的函数
来拟合这组数据,要求所得的拟合曲线能最好的反映数据的基本趋势(即使
最好地逼近
,而不必满足插值原则。
因此没必要取
=
,只要使
尽可能地小)。
原理:
给定数据点
。
求近似曲线
。
并且使得近似曲线与
的偏差最小。
近似曲线在该点处的偏差
,i=1,2,...,m。
常见的曲线拟合方法:
1.使偏差绝对值之和最小
2.使偏差绝对值最大的最小
3.使偏差平方和最小
最小二乘法:
按偏差平方和最小的原则选取拟合曲线,并且采取二项式方程为拟合曲线的方法,称为最小二乘法。
推导过程:
1.设拟合多项式为:
2.各点到这条曲线的距离之和,即偏差平方和如下:
3.问题转化为求待定系数
...
对等式右边求
偏导数,因而我们得到了:
.......
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.15560
图1最小二乘法曲线拟合示例
对比检验拟合的有效性:
例:
在[0,π]区间上对正弦函数进行拟合,然后在[0,2π]区间画出图形,比较拟合区间和非拟合区间的图形,考察拟合的有效性。
在MATLAB中输入如下代码:
clear
x=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.0000
mu=
R:
[10x10double]
df:
22
normr:
1.6178e-07
MATLAB的最优化工具箱还提供了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.8955
0.10003.5639
0.20002.5173
0.30001.9790
0.40001.8990
0.50001.3938
0.60001.1359
0.70001.0096
0.80001.0343
0.90000.8435
1.00000.6856
1.10000.6100
1.20000.5392
1.30000.3946
1.40000.3903
1.50000.5474
1.60000.3459
1.70000.1370
1.80000.2211
1.90000.1704
2.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)
holdon
plot(t,F(x,t))
holdoff
Fsumsquares=@(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)
typefitvector
x02=[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)
holdon
plot(t,F(xbad,t),'g')
legend('Data','Globalfit','Badlocalfit','Location','NE')
holdoff
fprintf(['Theresidualnormatthegoodendingpointis%f,'...
'andtheresidualnormatthebadendingpointis%f.'],...
resnorm,resnormbad)
displayEndOfDemoMessage(m)
拟合效果如下:
直线的最小二乘拟合:
y=a+bx
式中有两个待定参数,a代表截距,b代表斜率。
对于等精度测量所得到的N组数据(xi,yi),i=1,2……,N,xi值被认为是准确的,所有的误差只联系着yi。
下面利用最小二乘法把观测数据拟合为直线。
用最小二乘法估计参数时,要求观测值yi的偏差的加权平方和为最小。
对于等精度观测值的直线拟合来说,可使下式的值最小:
上式分别对a、b求偏导得:
整理后得到方程组:
解上述方程组便可求得直线参数a和b的最佳估计值。
1、可看成是一阶多项式拟合,跟前面曲线拟合方法一样。
2、利用linefit()函数进行最小二乘的直线拟合
使用:
clear
x=[0.511.522.53];
y=[1.752.453.814.888.6];
[k,b]=linefit(x,y)%得到斜率k和常数b
y1=polyval([k,b],x);
plot(x,y1,’k-’,x,y,’k*’)
MATLAB一元到多元线性回归方程的计算和检验_XX文库
(研究生数理统计)多元线性回归及显著性检验Matlab程序(完美版)_XX文库