《matlab最小二乘法拟合.doc》由会员分享，可在线阅读，更多相关《matlab最小二乘法拟合.doc（2页珍藏版）》请在三一文库上搜索。
1、最小二乘法在曲线拟合中比较普遍。 拟合的模型主要有 1. 直线型 2. 多项式型 3. 分数函数型 4. 指数函数型 5. 对数线性型 6. 高斯函数型 一般对于 LS问题，通常利用反斜杠运算“”、fminsearch 或优化工具箱提供的极小化函数求解。 在 Matlab 中，曲线拟合工具箱也提供了曲线拟合的图形界面操作。 在命 令提示符后键入： cftool ，即可根据数据， 选 择适当的拟合模型。 “ ”命令 命 令1.假设要拟合的多项式是：y=a+b*x+c*x2.首先建立设计矩阵X:X=o nes(size(x) x xA2;执行：para=Xy para中包含了三个参数：para(1
2、)=a;para(2)=b;para(3)=c; 这种方法对于系数是线性的模型也适应。 2. 假设要拟合： y=a+b*exp(x)+cx*exp(xA2) 设计矩阵 X 为 X=ones(size(x) exp(x) x.*exp(x.A2); para=Xy 3. 多重回归(乘积回归) 设要拟合： y=a+b*x+c*t ，其中 x 和 t 是预测变量， y 是响应变量。设计矩阵为 X=ones(size(x) x t %注意 x,t大小相等！ para=Xy polyfit 函数 polyfit函数不需要输入设计矩阵， 在参数估计中， polyfit 会根据输入的数据生成设计 矩阵。 1
3、. 假设要拟合的多项式是： y=a+b*x+c*xA2 p=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(
4、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
5、);% 初值 para=fminsearch(p) (y-f(p,x).A2,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
6、,x) 2.350176*p(1)*(1-1/p(2)*(1-(1-x(:,1).A(1/p(2).Ap(2).A2.*(x(:,1).A (-1/p(2)-1).A(-p(2).*x(:,1).A(-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 年的总人口情 况
7、的曲线拟合clear all;close all; %cftool 提供了可视化的曲线拟合！ t=1900 1910 19201930 1940 1950 1960 1970 1980 1990 2000' y=75.995 91.972 105.711123.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=
8、a+bx mx=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 2
9、020;% 预测年份 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
10、-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)42 %RA2 %拟合误差估计 CI YearAdd=t;Pre
11、Yea门；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
12、 plot(s,Y(:,col)-2*delta(:,col),'r:');plot(tPreYear,PrePo(:,col),'>');%预测值plot(tPreYear,PrePo(:,col)+2*Predelta(:,col);%预测 95% CIplot(tPreYear,PrePo(:,col)-2*Predelta(:,col);axis(-1.21.8 0 400);set(gca,'XTick',tYearAdd,'XTickLabel',sprintf('%d|',YearAdd);ti
13、tle(CFtitcol); hold off; end figure;%残差图 for col=1:nsubplot(2,2,col); plot(Y(:,i),res(:,i),'o'); end一个非线性的应用例子 (多元情况) 一个非线性的 在百度知道中，要拟合 y=a*x1An1+b*x2An2+c*x3An3 %注：只是作为应用，模型不一定正确！ %x2=x3!y=1080.94 1083.03 1162.80 1155.61 1092.82 1099.26 1161.06 1258.051299.03 1298.30 1440.22 1641.30 1672.21
14、 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.0251.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.1751.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' x=x1 x2 x3; f=(p,x)p(1)*x(:,1).Ap(2)+p(3)*x(:,2).Ap(4)+p(5)*x(:,3).Ap(6); p0=ones(6,1);失败p=fmi nsearch(p)sum(y-f(p,x)42,p0) res=y-f(p,x); res2=res.A2 % 的模型