matlab最小二乘法的非线性参数拟合

1、matlab最小二乘法的非线性参数拟合首先说一下匿名函数：在创建匿名函数时，Matlab记录了关于函数的信息,当使用句柄调用该函数的时候，Matlab不再进行搜索，而是立即执行该函数， 极大提高了效率。所以首选匿名函数。具体拟合时可以使用的方法如下： 1曲线拟合工具箱提供了很多拟合函数，使用简单 非线性拟合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.4

2、4 0.466 0.479 0.309'f=(p,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最小二乘法在曲线拟合中比较普遍。拟合的模型主要有1. 直线型2. 多项式型3. 分数函数型4. 指数函数

3、型5. 对数线性型6. 高斯函数型一般对于LS问题，通常利用反斜杠运算"”、fminsearch或优化工具箱提供的极小化函数求解。在 Matlab中，曲线拟合工具箱也提供了曲线拟合的图 形界面操作。在命令提示符后键入：cftool ,即可根据数据，选择适当的拟合模 型。“ ”命令1. 假设要拟合的多项式是：y=a+b*x+c*xA2.首先建立设计矩阵X： X=ones(size(x) x xA2;执行：para=Xypara 中包含了三个参数：para(1)=a;para(2)=b;para(3)=c;这种方法对于系数是线性的模型也适应。2. 假设要拟合：y=a+b*exp(x)+c

4、x*exp(xA2)设计矩阵X为X=ones(size(x) exp(x) x.*exp(x.A2);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*xA2p=polyfit(x,y,2)然后可以使用polyval在t处预测：y_hat=polyval(p,t)polyfit函数可

5、以给出置信区问。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函数进 行曲线拟合。假设要

6、拟合：y = a+b*exp(x)+c*exp(x.?2)首先建立函数，可以通过 m文件或函数句柄建立：x=.'；y=.'f=(p,x) p(1)+p (2) *exp(x)+p *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).A2,p0) % 可以输出 Hessian 矩阵res=y-f(para,x)% 拟合残差3. 多项式型的一个例子1900-2000年的总人口情况的曲线拟合clear all;

7、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=

8、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=zero

9、s(length(t),n);res=zeros(size(Y);delta=zeros(size(Y);PrePo=zeros(length(PreYear),n);Predelta=zeros(size(PrePo);for i=1:np S(i)=polyfit(s,y,i);Y(:,i) delta(:,i)=polyval(p,s,S(i);%拟合的 YPrePo(:,i) Predelta(:,i)=polyval(p,tPreYear,S(i);%预测res(:,i)=y-Y(:,i);%残差end% plot(s,Y);%2009a自动添加不同颜色%legend('da

10、ta','regression line','1st poly','2nd poly','3rd poly','4thpoly',2)% plot(tPreYear,PrePo,'>');% hold off% plot(Y,res,'o');% 残差图r=corr(s,Y).A2 %RA2以拟合误差估计CIYearAdd=t;PreYear'tYearAdd=s;tPreYear'CFtit=' 一阶拟合','二阶拟合'

11、,'三阶拟合','四阶拟合'for col=1:nsubplot(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% CIplot(s,Y(:,col)-2*delta(:,col),'r:');plot(tPreYear,PrePo(:,col),'

12、>');%预测值plot(tPreYear,PrePo(:,col)+2*Predelta(:,col);%预测 95% CIplot(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:nsubplot(2,2,col);plot(Y(:,i),

13、res(:,i),'o');end4非线性的应用例子(多元情况)要拟合 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 1612.73 1658.64 1752.42 1837.992099.29 2675.47 2786.33 2881.07'x1=1 1.05 1.1 1.15 1.2 1.25 1.3

14、 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.325 1.350 1.375 1.4 1.425 1.45 1.475 1.5'x=

15、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=fminsearch(p)sum(y-f(p,x).A2,p0)res=y-f(p,x);res2=res.A2 %失败的模型Matlab自定义函数自定义函数的途径：M文件函数(M file function)在线函数(Inline Function)匿名函数(Anonymous Function)1.M文件函数范例function c=myadd(a,b)以这里可以写函数的使用说明，前面以 册头以拄工作区中，help

16、 myadd将显示此处的说明c=a+b;%end %必须的第一行function 告诉Matlab这是一个函数,a,b是输入,c是输出,myadd是 函数名。以m文件定义的函数必须保存为函数名的形式， 上例中，函数应保存为 myadd.m要使用myadd函数，该函数必须在 Matlab的搜索路径中。调用方式：在Matlab命令符后输入a=1;b=2;c=myadd(a,b)2. 在线函数(Inline Function)通常作为参数传递给另外一个函数。比如fminsearch , lsqcurvefit 等函数需要以函数作为参数。在线函数从字符申表达式创建函数，例如：f=inline('x.A2,x,);创建了函数f(x)=xA2。要计算f(3)，在工作区输入f(3)即可。f(2 3 4) 计算 在x=2 3 4时的值 f=inline('x+y','x','y')创建了二元函数f(x,y)=x+y ，工作区输入f(2,3)计算2+3,等同于f(f,2,3)。3. 匿名函数(Anonymous Function)匿名函数使用函数句柄来表示匿名函数，定义形式为函数句柄=奁量名)函数表达式例如：f=(x) x.A2定义了函数f(x)=xA2 , f(2)计算在x=2处的值。匿名函数可以调