曲线拟合的Matlab实现

曲线拟合算法之MATLAB实现functionretstr=FitDrawpic(ModelNo,NetPara,SeqData,TargetData,DataDir)NNTWARNOFF%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ModelNo='1';%DataDir='.'; %表示当前目录%NetPara(1)=11;%NetPara(2)=3;%NetPara(3)=50;%NetPara(4)=1900;%NetPara(5)=2020;%NetPara(6)=70;%NetPara(7)=350;%NetPara(8)=2010;%SeqData=[1900.0000000000000001910.0000000000000001920.0000000000000001930.0000000000000001940.0000000000000001950.0000000000000001960.0000000000000001970.0000000000000001980.0000000000000001990.0000000000000002000.000000000000000];%TargetData=[75.9950000091.97200000105.71100000123.20300000131.66900000150.69700000179.32300000203.21200000226.50500000249.63300000281.42200000];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%保留原目录olddir=pwd;%进入数据所在目录cd(DataDir);InputDim=NetPara(1);cd(DataDir);InputDim=NetPara(1);FitFact=NetPara(2);DataNum=NetPara(3);%样本组数%拟合阶数%拟合数据点数%%外推最小值%外推最大值MinSeqValue=NetPara(4);MaxSeqValue=NetPara(5);MinValue=NetPara(6);%输出最小值MaxValue=NetPara(7);%输出最大值FitPoint=NetPara(8); %预测点x=SeqData;y=TargetData;newx=linspace(MinSeqValue,MaxSeqValue,DataNum);%读取拟合曲线系数fcoef=fopen(sprintf('coef%s%s',ModelNo,'.dat'),'r');[coef,count]=fscanf(fcoef,'%f',[1,inf]);fclose(fcoef);z=polyval(coef,newx)';pz=polyval(coef,FitPoint);holdonaxis([MinSeqValue,MaxSeqValue,MinValue,MaxValue]);plot(newx,z,x,y,'o');plot(FitPoint,pz,'ks');text(FitPoint,pz+15,num2str(pz));holdofftitle]曲线拟合示意图,阶数=',int2str(FitFact)])xlabel('序列值'),ylabel('输出值'),gridcd(olddir);retstr=pz;functionretstr=FitForecast(ModelNo,InputData,DataDir)NNTWARNOFF%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ModelNo='1';%DataDir='.'; %表示当前目录%InputData=2010;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%保留原目录olddir=pwd;%进入数据所在目录cd(DataDir);%读取拟合曲线系数fcoef=fopen(sprintf('coef%s%s',ModelNo,'.dat'),'r');[coef,count]=fscanf(fcoef,'%f',[1,inf]);fclose(fcoef);pz=polyval(coef,InputData);cd(olddir)retstr=pz;

functionretstr=FitPoly(ModelNo,NetPara,SeqData,TargetData,DataDir)NNTWARNOFF%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ModelNo='1';%DataDir='.'; %表示当前目录%NetPara(1)=11;%NetPara(2)=3;%NetPara(3)=100;%SeqData=(1900:10:2000);%TargetData=[75.99591.972105.711123.203131.669...%150.697179.323203.212226.505249.633281.422];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%保留原目录olddir=pwd;%进入数据所在目录cd(DataDir);InputDim=NetPara(1);FitFact=NetPara(2);DataNum=NetPara(3);x=SeqData;y=TargetData;InputDim=NetPara(1);FitFact=NetPara(2);DataNum=NetPara(3);x=SeqData;y=TargetData;%样本组数%拟合阶数%拟合据数点数%[P,S,MU]=POLYFIT(X,Y,N)findsthecoefficientsofapolynomial%inXHAT=(X-MU(1))/MU(2)whereMU(1)=mean(X)andMU(2)=std(X).%Thiscenteringandscalingtransformationimprovesthenumerical%propertiesofboththepolynomialandthefittingalgorithm.%ThestructureScontainstheCholeskyfactoroftheVandermonde%matrix(R),thedegreesoffreedom(df),andthenormofthe%residuals(normr)asfields.[coef,s,mu]=polyfit(x,y,FitFact);[coef]=polyfit(x,y,FitFact);%将参数值写入文件coefficientfcoef=fopen(sprintf('coef%s%s',ModelNo,'.dat'),'w');fsnormr=fopen(sprintf('snormr%s%s',ModelNo,'.dat'),'w');fxhat=fopen(sprintf('xhat%s%s',ModelNo,'.dat'),'w');%多项式拟合系数fprintf(fcoef,'%25.15f',coef);%残差fprintf(fsnormr,'%15.8f',s.normr);%各点误差fprintf(fxhat,'%10.6f',(x-mu(1))/mu(2));fclose(fcoef);fclose(fsnormr);fclose(fxhat);cd(olddir);retstr=1;functionretstr=FitSimu(ModelNo,NetPara,SeqData,TargetData,DataDir)NNTWARNOFF%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ModelNo='1';%DataDir='.'; %表示当前目录%NetPara(1)=11;%NetPara(2)=3;%NetPara(3)=50;%NetPara(4)=1900;%NetPara(5)=2020;%NetPara(6)=70;%NetPara(7)=350;%NetPara(8)=2010;%SeqData=[1900.0000000000000001910.0000000000000001920.0000000000000001930.0000000000000001940.0000000000000001950.0000000000000001960.0000000000000001970.0000000000000001980.0000000000000001990.0000000000000002000.000000000000000];%TargetData=[75.9950000091.97200000105.71100000123.20300000131.66900000150.69700000179.32300000203.21200000226.50500000249.63300000281.42200000];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%保留原目录olddir=pwd;%进入数据所在目录cd(DataDir);InputDim=NetPara(1); %样本组数FitFact=NetPara(2); %拟合阶数DataNum=NetPara(3); %拟合数据点数MinSeqValue=NetPara(4);%外推最小值MaxSeqValue=NetPara(5);%外推最大值MinValue=NetPara(6);%输出最小值MaxValue=NetPara(7);%输出最大值FitPoint=NetPara(8); %预测点%读取拟合曲线系数fcoef=fopen(sprintf('coef%s%s',ModelNo,'.dat'),'r');[coef,count]=fscanf(fcoef,'%f',[1,inf]);fclose(fcoef);pz=polyval(coef,FitPoint);cd(olddir);retstr=pz;曲线拟合与插值在大量的应用领域中，人们经常面临用一个解析函数描述数据(通常是测量值)的任务。对这个问题有两种方法。在插值法里，数据假定是正确的，要求以某种方法描述数据点之间所发生的情况。这种方法在下一节讨论。这里讨论的方法是曲线拟合或回归。人们设法找出某条光滑曲线，它最佳地拟合数据，但不必要经过任何数据点。图11.1说明了这两种方法。标有'o'的是数据点；连接数据点的实线描绘了线性内插，虚线是数据的最佳拟合。11.1曲线拟合曲线拟合涉及回答两个基本问题：最佳拟合意味着什么？应该用什么样的曲线？可用许多不同的方法定义最佳拟合，并存在无穷数目的曲线。所以，从这里开始，我们走向何方？正如它证实的那样，当最佳拟合被解释为在数据点的最小误差平方和，且所用的曲线限定为多项式时，那么曲线拟合是相当简捷的。数学上，称为多项式的最小二乘曲线拟合。如果这种描述使你混淆，再研究图11.1。虚线和标志的数据点之间的垂直距离是在该点的误差。对各数据点距离求平方，并把平方距离全加起来，就是误差平方和。这条虚线是使误差平方和尽可能小的曲线，即是最佳拟合。最小二乘这个术语仅仅是使误差平方和最小的省略说法。在MATLAB中，函数polyfit求解最小二乘曲线拟合问题。为了阐述这个函数的用法，让我们以上面图11.1中的数据开始。»x=[0.1.2.3.4.5.6.7.8.91];»y=[-.4471.9783.286.167.087.347.669.569.489.3011.2];为了用polyfit，我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶次或度。如果我们选择n=1作为阶次，得到最简单的线性近似。通常称为线性回归。相反，如果我们选择n=2作为阶次，得到一个2阶多项式。现在，我们选择一个2阶多项式。»n=2;%polynomialorder»p=polyfit(x,y,n)p=-9.810820.1293-0.0317polyfit的输出是一个多项式系数的行向量。其解是y=—9.8108x2+20.1293x—0.0317。为了将曲线拟合解与数据点比较，让我们把二者都绘成图。»xi=linspace(0,1,100);%x-axisdataforplotting»z=polyval(p,xi);为了计算在xi数据点的多项式值，调用MATLAB的函数polyval。»plot(x,y,'o',x,y,xi,z,':')画出了原始数据x和y,用'o'标出该数据点，在数据点之间，再用直线重画原始数据，并用点'：'线，画出多项式数据xi和z。»xlabel('x'),ylabel('y=f(x)'),title('SecondOrderCurveFitting')将图作标志。这些步骤的结果表示于前面的图11.1中。多项式阶次的选择是有点任意的。两点决定一直线或一阶多项式。三点决定一个平方或2阶多项式。按此进行，n+1数据点唯一地确定n阶多项式。于是，在上面的情况下，有11个数据点，我们可选一个高达10阶的多项式。然而，高阶多项式给出很差的数值特性，人们不应选择比所需的阶次高的多项式。此外，随着多项式阶次的提高，近似变得不够光滑，因为较高阶次多项式在变零前，可多次求导。例如，选一个10阶多项式»pp=polyfit(x,y,10);»formatshorte%changedisplayformat»pp.'%displaypolynomialcoefficientsasacolumnans=-4.6436e+0052.2965e+006-4.8773e+0065.8233e+006-4.2948e+0062.0211e+006-6.0322e+0051.0896e+005-1.0626e+0044.3599e+002-4.4700e-001要注意在现在情况下，多项式系数的规模与前面的2阶拟合的比较。还要注意在最小(-4.4700e-001)和最大(5.8233e+006)系数之间有7个数量级的幅度差。将这个解作图，并把此图与原始数据及2阶曲线拟合相比较，结果如何呢？»zz=polyval(pp,xi);%evaluate10thorderpolynomial»plot(x,y,'o',xi,z,'：',xi,zz)%plotdata»xlabel('x'),ylabel('y=f(x)'),title('2ndand10thOrdercurveFitting')在下面的图11.2中，原始数据标以'o',2阶曲线拟合是虚线，10阶拟合是实线。注意，在10阶拟合中，在左边和右边的极值处，数据点之间出现大的纹波。当企图进行高阶曲线拟合时，这种纹波现象经常发生。根据图11.2，显然，‘越多就越好'的观念在这里不适用。11.2一维插值正如在前一节对曲线拟合所描述的那样，插值定义为对数据点之间函数的估值方法，这些数据点是由某些集合给定。当人们不能很快地求出所需中间点的函数值时，插值是一个有价值的工具。例如，当数据点是某些实验测量的结果或是过长的计算过程时，就有这种情况。或许最简单插值的例子是MATLAB的作图。按缺省，MATLAB用直线连接所用的数据点以作图。这个线性插值猜测中间值落在数据点之间的直线上。当然，当数据点个数的增加和它们之间距离的减小时，线性插值就更精确。例如，»x1=linspace(0,2*pi,60);»x2=linspace(0,2*pi,6);»plot(x1,sin(x1),x2,sin(x2),'-')»xlabel('x'),ylabel('sin(x)'),title('LinearInterpolation')如曲线拟合一样，插值要作决策。根据所作的假设，有多种插值。而且，可以在一维以上空间中进行插值。即如果有反映两个变量函数的插值，z=f(x,y),那么就可在x之间和在y之间，找出z的中间值进行插值。MATLAB在一维函数interpl和在二维函数interp2中，提供了许多的插值选择。其中的每个函数将在下面阐述。为了说明一维插值，考虑下列问题，12小时内，一小时测量一次室外温度。数据存储在两个MATLAB变量中。»hours=1:12;%indexforhourdatawasrecorded»temps=[589152529313022252724];%recordedtemperatures»plot(hours,temps,hours,temps,'+')%viewtemperatures»title('Temperature')»xlabel('Hour'),ylabel('DegreesCelsius')MATLAB画出了数据点线性插值的直线。为了计算在任意给定时间的温度，人们可试着对可视的图作解释。另外一种方法，可用函数interpl。»t=interp1(hours,temps,9.3)%estimatetemperatureathour=9.3t=22.9000»t=interpl(hours,temps,4.7)%estimatetemperatureathour=4.7t=22»t=interpl(hours,temps,[3.26.57.lll.7])%findtempatmanypoints!t=l0.200030.000030.900024.9000interpl的缺省用法是由interp1(x,y,xo)来描述，这里x是独立变量(横坐标)，y是应变量(纵坐标)，xo是进行插值的一个数值数组。另外，该缺省的使用假定为线性插值。若不采用直线连接数据点，我们可采用某些更光滑的曲线来拟合数据点。最常用的方法是用一个3阶多项