matlab最小二乘法的非線(xiàn)性參數(shù)擬合

1、matlab 最小二乘法的非線(xiàn)性參數(shù)擬合首先說(shuō)一下匿名函數(shù)：在創(chuàng)建匿名函數(shù)時(shí)， Matlab 記錄了關(guān)于函數(shù)的信息， 當(dāng)使用句柄調(diào)用該函數(shù)的時(shí)候， Matlab 不再進(jìn)行搜索，而是立即執(zhí)行該函數(shù)， 極大提高了效率。所以首選匿名函數(shù)。具體擬合時(shí)可以使用的方法如下： 1 曲線(xiàn)擬合工具箱 提供了很多擬合函數(shù)，使用簡(jiǎn)單 非線(xiàn)性擬合 nlinfit 函數(shù) clear all; x1='x2='x=x1 x2; y= ' f=(p,x)*p(1)*(1-1/p (2) )*(1-(1-x(:,1)./p (2)邛(2)A2.*(x(:,1)A(-1/p (2)-1)A(-p(2).

2、*x(:,1).(1/p (2).*x(:,2);p0=8 'opt=optimset('TolFun',1e-3,'TolX',1e-3);%p R=nlinfit(x,y,f,p0,opt)2 最小二乘法在曲線(xiàn)擬合中比較普遍。 擬合的模型主要有1. 直線(xiàn)型2. 多項(xiàng)式型3. 分?jǐn)?shù)函數(shù)型4. 指數(shù)函數(shù)型5. 對(duì)數(shù)線(xiàn)性型6. 高斯函數(shù)型一般對(duì)于LS問(wèn)題，通常利用反斜杠運(yùn)算 “”、fminsearch或優(yōu)化工具箱提供 的極小化函數(shù)求解。在 Matlab 中，曲線(xiàn)擬合工具箱也提供了曲線(xiàn)擬合的圖形界 面操作。在命令提示符后鍵入： cftool ，即可根據(jù)數(shù)據(jù)，

3、選擇適當(dāng)?shù)臄M合模型?！啊泵?. 假設(shè)要擬合的多項(xiàng)式是：y=a+b*x+c*x9.首先建立設(shè)計(jì)矩陣X：X=o nes(size(x) x xA2;執(zhí)行：para=Xypara 中包含了三個(gè)參數(shù)： para(1)=a;para(2)=b;para(3)=c; 這種方法對(duì)于系數(shù)是線(xiàn)性的模型也適應(yīng)。2. 假設(shè)要擬合： y=a+b*exp(x)+cx*exp(xA2)設(shè)計(jì)矩陣 X 為X=ones(size(x) exp(x) x.*exp(x.A2);para=Xy3. 多重回歸 (乘積回歸 )設(shè)要擬合：y=a+b*x+c*t，其中x和t是預(yù)測(cè)變量，y是響應(yīng)變量。設(shè)計(jì)矩陣為X=ones(size(x)

4、 x t %注意 x,t 大小相等！para=Xypolyfit 函數(shù)polyfit 函數(shù)不需要輸入設(shè)計(jì)矩陣，在參數(shù)估計(jì)中， polyfit 會(huì)根據(jù)輸入的數(shù)據(jù)生成 設(shè)計(jì)矩陣。1. 假設(shè)要擬合的多項(xiàng)式是：y=a+b*x+c*xA2p=polyfit(x,y,2)然后可以使用 polyval 在 t 處預(yù)測(cè)：y_hat=polyval(p,t)polyfit 函數(shù)可以給出置信區(qū)間。p S=polyfit(x,y,2) %5中包含了標(biāo)準(zhǔn)差y_fit,delta = polyval(p,t,S) %按照擬合模型在 t 處預(yù)測(cè)在每個(gè) t 處的 95%CI為:*delta, y_fit+*delta)2.

5、 指數(shù)模型也適應(yīng)假設(shè)要擬合： y = a+b*exp(x)+c*exp(x.2)p=polyfit(x,log(y),2)fminsearch 函數(shù)fmin search是優(yōu)化工具箱的極小化函數(shù)。LS問(wèn)題的基本思想就是殘差的平方和(一 種范數(shù)，由此，LS產(chǎn)生了許多應(yīng)用)最小，因此可以利用fminsearch函數(shù)進(jìn)行曲 線(xiàn)擬合。假設(shè)要擬合： y = a+b*exp(x)+c*exp(x.2) 首先建立函數(shù)，可以通過(guò) m 文件或函數(shù)句柄建立： x='y='f=(p,x) p(1)+p(2)*exp(x)+p(3)*exp(x.2) %注意向量化 :p(1)=a;p(2)=b;p(3

6、)=c; %可以根據(jù)需要選擇是否優(yōu)化參數(shù)%opt=options()pO=o nes(3,1);%初值para=fminsearch(p) (y-f(p,x).A2,pO) %可以輸出 Hessian矩陣 res=y-f(para,x)%擬合殘差3. 多項(xiàng)式型 的一個(gè)例子1900-2000年的總?cè)丝谇闆r的曲線(xiàn)擬合clear all;close all;%cftool 提供了可視化的曲線(xiàn)擬合！t=1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000'y='%t太大，以t的幕作為基函數(shù)會(huì)導(dǎo)致設(shè)計(jì)矩陣尺度太差，列變量幾乎線(xiàn)性相依

7、。 變換為 -1 1上s=(t-1950)/50;%plot(s,y,'ro'); %回歸線(xiàn)： y=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(' 多項(xiàng)式擬合 '); set(gca,'XTick',s,'XTickLabel',spri

8、ntf('%d|',t);%hold on;n=4;PreYear=2010 2015 2020;%預(yù)測(cè)年份 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:np S(i)=polyfit(s,y,i);Y(:,i) delta(:,i)=polyval(p,s,S(i);%合的 Y PrePo(:,i) Predelta

9、(:,i)=polyval(p,tPreYear,S(i);預(yù)測(cè)res(:,i)=y-Y(:,i);%殘差end% plot(s,Y);%2009a自動(dòng)添加不同顏色% 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)

10、.A2 %RA2%擬合誤差估計(jì) CIYearAdd=t;PreYear'tYearAdd=s;tPreYear'CFtit=一階擬合','二階擬合','三階擬合','四階擬合'for col=1:nsubplot(3,2,col+2);plot(s,y,'ro',s,Y(:,col),'g-');%原 始數(shù)據(jù)和擬合數(shù)據(jù)legend('Original','Fitted',2);hold on;plot(s,Y(:,col)+2*delta(:,col),

11、9;r:');%95% CI plot(s,Y(:,col)-2*delta(:,col),'r:');plot(tPreYear,PrePo(:,col),'>');%預(yù)測(cè)值plot(tPreYear,PrePo(:,col)+2*Predelta(:,col);%預(yù)測(cè) 95% CI plot(tPreYear,PrePo(:,col)-2*Predelta(:,col);axis( 0 400); set(gca,'XTick',tYearAdd,'XTickLabel',sprintf('%d|'

12、;,YearAdd); title(CFtitcol);hold off;endfigure;%殘差圖for col=1:nsubplot(2,2,col);plot(Y(:,i),res(:,i),'o');end4 非線(xiàn)性的 應(yīng)用例子 ( 多元情況 )要擬合 y=a*xMn1+b*x2"n2+c*x3"n3%注：只是作為應(yīng)用，模型不一定正確！%x2=x3!y='x1=12'x2=1'x3=1'x=x1 x2 x3;f=(p,x) p(1)*x(:,1)Ap(2)+p(3)*x(:,2)Ap(4)+p(5)*x(:,3)Ap

13、(6);p0=ones(6,1);p=fmi nsearch(p)sum(y-f(p,x)42,pO)res=y-f(p,x);res2=res.A2 %失敗的模型Matlab 自定義函數(shù)自定義函數(shù)的途徑：M 文件函數(shù) (M file function)在線(xiàn)函數(shù) (Inline Function)匿名函數(shù) (Anonymous Function)文件函數(shù)范例function c=myadd(a,b)%這里可以寫(xiě)函數(shù)的使用說(shuō)明，前面以 %開(kāi)頭 %在工作區(qū)中， help myadd 將顯示此處的說(shuō)明 c=a+b;%end %非必須的第一行 function 告訴 Matlab 這是一個(gè)函數(shù)， a,

14、b 是輸入， c 是輸出， myadd 是函 數(shù)名。以 m 文件定義的函數(shù)必須保存為函數(shù)名的形式，上例中，函數(shù)應(yīng)保存為。 要使用 myadd 函數(shù)，該函數(shù)必須在 Matlab 的搜索路徑中。調(diào)用方式：在 Matlab 命令符后輸入 a=1;b=2;c=myadd(a,b)2.在線(xiàn)函數(shù) (Inline Function)通常作為參數(shù)傳遞給另外一個(gè)函數(shù)。比如fminsearch， lsqcurvefit 等函數(shù)需要以函數(shù)作為參數(shù)。在線(xiàn)函數(shù)從字符串表達(dá)式創(chuàng)建函數(shù)，例如：f=in li ne('x.A2','x');創(chuàng)建了函數(shù)f(x)=xA20要計(jì)算f(3)，在工作區(qū)輸入f(3)即可。f(2 3 4)計(jì)算在x=2 34 時(shí)的值f=inline('x+y','x','y')創(chuàng)建了二元函數(shù)f(x,y)=x+y,工作區(qū)輸入f(2,3)計(jì)算2+3,等同于f(f,2,3)。3.匿名函數(shù) (Anonymous Function) 匿名函數(shù)使用函數(shù)句柄來(lái)表示匿名函數(shù)，定義形式為 函數(shù)句柄 =(變量名) 函數(shù)表達(dá)式例如：f=(x) x.A2定義了函數(shù)f(x)=xA2, f(2)計(jì)算在x=2處的值。匿名函數(shù)可以調(diào)用 Matlab 函數(shù)，也可以使用工作區(qū)中存在的變量，例如 a=2;f=(x) x.A