Matlab最小二乘法曲线拟合.doc

最小二乘法在曲线拟合中比较普遍。拟合的模型主要有1.直线型2.多项式型3.分数函数型4.指数函数型5.对数线性型6.高斯函数型.一般对于LS问题，通常利用反斜杠运算“”、fminsearch或优化工具箱提供的极小化函数求解。在Matlab中，曲线拟合工具箱也提供了曲线拟合的图形界面操作。在命令提示符后键入：cftool，即可根据数据，选择适当的拟合模型。“”命令1.假设要拟合的多项式是：y=a+b*x+c*x2.首先建立设计矩阵X：X=ones(size(x) x x2;执行：para=Xypara中包含了三个参数：para(1)=a;para(2)=b;para(3)=c;这种方法对于系数是线性的模型也适应。2.假设要拟合：y=a+b*exp(x)+cx*exp(x2)设计矩阵X为X=ones(size(x) exp(x) x.*exp(x.2);para=Xy3.多重回归(乘积回归)设要拟合：y=a+b*x+c*t，其中x和t是预测变量，y是响应变量。设计矩阵为X=ones(size(x) x t %注意x,t大小相等！para=Xypolyfit函数polyfit函数不需要输入设计矩阵，在参数估计中，polyfit会根据输入的数据生成设计矩阵。1.假设要拟合的多项式是：y=a+b*x+c*x2p=polyfit(x,y,2)然后可以使用polyval在t处预测：y_hat=polyval(p,t)polyfit函数可以给出置信区间。p S=polyfit(x,y,2) %S中包含了标准差y_fit,delta = polyval(p,t,S) %按照拟合模型在t处预测在每个t处的95%CI为：(y_fit-1.96*delta, y_fit+1.96*delta)2.指数模型也适应假设要拟合：y = a+b*exp(x)+c*exp(x.?2)p=polyfit(x,log(y),2)fminsearch函数fminsearch是优化工具箱的极小化函数。LS问题的基本思想就是残差的平方和(一种范数，由此，LS产生了许多应用)最小，因此可以利用fminsearch函数进行曲线拟合。假设要拟合：y = a+b*exp(x)+c*exp(x.?2)首先建立函数，可以通过m文件或函数句柄建立：x=.;y=.;f=(p,x) p(1)+p(2)*exp(x)+p(3)*exp(x.?2) %注意向量化:p(1)=a;p(2)=b;p(3)=c;%可以根据需要选择是否优化参数%opt=options()p0=ones(3,1);%初值para=fminsearch(p) (y-f(p,x).2,p0) %可以输出Hessian矩阵res=y-f(para,x)%拟合残差曲线拟合工具箱提供了很多拟合函数，对大样本场合比较有效！非线性拟合nlinfit函数clear all;x1=0.4292 0.4269 0.381 0.4015 0.4117 0.3017;x2=0.00014 0.00059 0.0126 0.0061 0.00425 0.0443;x=x1 x2;y=0.517 0.509 0.44 0.466 0.479 0.309;f=(p,x) 2.350176*p(1)*(1-1/p(2)*(1-(1-x(:,1).(1/p(2).p(2).2.*(x(:,1).(-1/p(2)-1).(-p(2).*x(:,1).(-1/p(2)-0.5).*x(:,2);p0=8 0.5;opt=optimset(TolFun,1e-3,TolX,1e-3);%p R=nlinfit(x,y,f,p0,opt) 例子例子例子例子例子例子例子例子例子例子例子例子例子例子例子例子直线型例子2.多项式型的一个例子1900-2000年的总人口情况的曲线拟合clear all;close all;%cftool提供了可视化的曲线拟合！t=1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000;y=75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 281.4220;%t太大，以t的幂作为基函数会导致设计矩阵尺度太差，列变量几乎线性相依。变换为-1 1上s=(t-1950)/50;%plot(s,y,ro);%回归线：y=a+bxmx=mean(s);my=mean(y);sx=std(s);sy=std(y);r=corr(s,y);b=r*sy/sx;a=my-b*mx;rline=a+b.*s;figure;subplot(3,2,1 2)plot(s,y,ro,s,rline,k);%title(多项式拟合);set(gca,XTick,s,XTickLabel,sprintf(%d|,t);%hold on;n=4;PreYear=2010 2015 2020;%预测年份tPreYear=(PreYear-1950)/50;Y=zeros(length(t),n);res=zeros(size(Y);delta=zeros(size(Y);PrePo=zeros(length(PreYear),n);Predelta=zeros(size(PrePo);for i=1:n p S(i)=polyfit(s,y,i); Y(:,i) delta(:,i)=polyval(p,s,S(i);%拟合的Y PrePo(:,i) Predelta(:,i)=polyval(p,tPreYear,S(i);%预测 res(:,i)=y-Y(:,i);%残差end% plot(s,Y);%2009a自动添加不同颜色% legend(data,regression line,1st poly,2nd poly,3rd poly,4th poly,2)% plot(tPreYear,PrePo,);% hold off% plot(Y,res,o);%残差图r=corr(s,Y).2 %R2%拟合误差估计CIYearAdd=t;PreYear;tYearAdd=s;tPreYear;CFtit=一阶拟合,二阶拟合,三阶拟合,四阶拟合;for col=1:n subplot(3,2,col+2); plot(s,y,ro,s,Y(:,col),g-);%原始数据和拟合数据 legend(Original,Fitted,2); hold on; plot(s,Y(:,col)+2*delta(:,col),r:);%95% CI plot(s,Y(:,col)-2*delta(:,col),r:); plot(tPreYear,PrePo(:,col),);%预测值 plot(tPreYear,PrePo(:,col)+2*Predelta(:,col);%预测95% CI plot(tPreYear,PrePo(:,col)-2*Predelta(:,col); axis(-1.2 1.8 0 400); set(gca,XTick,tYearAdd,XTickLabel,sprintf(%d|,YearAdd); title(CFtitcol); hold off;endfigure;%残差图for col=1:n subplot(2,2,col); plot(Y(:,i),res(:,i),o);end一个非线性的应用例子(多元情况)在百度知道中，要拟合y=a*x1n1+b*x2n2+c*x3n3%注：只是作为应用，模型不一定正确！%x2=x3!y=1080.94 1083.03 1162.80 1155.61 1092.82 1099.26 1161.06 1258.05 1299.03 1298.30 1440.22 1641.30 1672.21 1612.73 1658.64 1752.42 1837.99 2099.29 2675.47 2786.33 2881.07; x1=1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5 1.55 1.6 1.65 1.7 1.75 1.8 1.85 1.9 1.95 2; x2=1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.250 1.275 1.3 1.325 1.350 1.375 1.4 1.425 1.45 1.475 1.5; x3=1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.250 1.275 1.3 1.3