matlab多元非線性回歸_第1頁
matlab多元非線性回歸_第2頁
matlab多元非線性回歸_第3頁
matlab多元非線性回歸_第4頁
matlab多元非線性回歸_第5頁
已閱讀5頁,還剩25頁未讀, 繼續(xù)免費閱讀

下載本文檔

版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請進行舉報或認(rèn)領(lǐng)

文檔簡介

matlab回歸(擬合)總結(jié)前言1、學(xué)三條命令polyfit(x,y,n)---擬合成一元冪函數(shù)(一元多次)regress(y,x)----可以多元,nlinfit(x,y,’fun’,beta0)(可用于任何類型的函數(shù),任意多元函數(shù),應(yīng)用范圍最主,最萬能的)2、同一個問題,這三條命令都可以使用,但結(jié)果肯定是不同的,因為擬合的近似結(jié)果,沒有唯一的標(biāo)準(zhǔn)的答案。相當(dāng)于咨詢多個專家。3、回歸的操作步驟:根據(jù)圖形(實際點),選配一條恰當(dāng)?shù)暮瘮?shù)形式(類型)---需要數(shù)學(xué)理論與基礎(chǔ)和經(jīng)驗。(并寫出該函數(shù)表達(dá)式的一般形式,含待定系數(shù))------選用某條回歸命令求出所有的待定系數(shù)。所以可以說,回歸就是求待定系數(shù)的過程(需確定函數(shù)的形式)一、回歸命令一元多次擬合polyfit(x,y,n);一元回歸polyfit;多元回歸regress---nlinfit(非線性)二、多元回歸分析對于多元線性回歸模型(其實可以是非線性,它通用性極高):設(shè)變量的n組觀測值為記,,則的估計值為排列方式與線性代數(shù)中的線性方程組相同(),擬合成多元函數(shù)---regress使用格式:左邊用b=[b,bint,r,rint,stats]右邊用=regress(y,x)或regress(y,x,alpha)---命令中是先y后x,---須構(gòu)造好矩陣x(x中的每列與目標(biāo)函數(shù)的一項對應(yīng))---并且x要在最前面額外添加全1列/對應(yīng)于常數(shù)項---y必須是列向量---結(jié)果是從常數(shù)項開始---與polyfit的不同。)其中:b為回歸系數(shù),的估計值(第一個為常數(shù)項),bint為回歸系數(shù)的區(qū)間估計,r:殘差,rint:殘差的置信區(qū)間,stats:用于檢驗回歸模型的統(tǒng)計量,有四個數(shù)值:相關(guān)系數(shù)r2、F值、與F對應(yīng)的概率p和殘差的方差(前兩個越大越好,后兩個越小越好),alpha:顯著性水平(缺省時為0.05,即置信水平為95%),(alpha不影響b,只影響bint(區(qū)間估計)。它越小,即置信度越高,則bint范圍越大。顯著水平越高,則區(qū)間就越?。ǚ祷匚鍌€結(jié)果)---如有n個自變量-有誤(n個待定系數(shù)),則b中就有n+1個系數(shù)(含常數(shù)項,---第一項為常數(shù)項)(b---b的范圍/置信區(qū)間---殘差r---r的置信區(qū)間rint-----點估計----區(qū)間估計此段上課時不要:----如果的置信區(qū)間(bint的第行)不包含0,則在顯著水平為時拒絕的假設(shè),認(rèn)為變量是顯著的.*******(而rint殘差的區(qū)間應(yīng)包含0則更好)。b,y等均為列向量,x為矩陣(表示了一組實際的數(shù)據(jù))必須在x第一列添加一個全1列。----對應(yīng)于常數(shù)項-------而nlinfit不能額外添加全1列。結(jié)果的系數(shù)就是與此矩陣相對應(yīng)的(常數(shù)項,x1,x2,……xn)。(結(jié)果與參數(shù)個數(shù):1/5=2/3-----y,x順序---x要額外添加全1列)而nlinfit:1/3=4------x,y順序---x不能額外添加全1列,---需編程序,用于模仿需擬合的函數(shù)的任意形式,一定兩個參數(shù),一為系數(shù)數(shù)組,二為自變量矩陣(每列為一個自變量)有n個變量---不準(zhǔn)確,x中就有n列,再添加一個全1列(相當(dāng)于常數(shù)項),就變?yōu)閚+1列,則結(jié)果中就有n+1個系數(shù)。x需要經(jīng)過加工,如添加全1列,可能還要添加其他需要的變換數(shù)據(jù)。相關(guān)系數(shù)r2越接近1,說明回歸方程越顯著;(r2越大越接近1越好)F越大,說明回歸方程越顯著;(F越大越好)與F對應(yīng)的概率p越小越好,一定要P<a時拒絕H0而接受H1,即回歸模型成立。乘余(殘差)標(biāo)準(zhǔn)差(RMSE)越小越好(此處是殘差的方差,還沒有開方)(前兩個越大越好,后兩個越小越好)regress多元(可通過變形而適用于任意函數(shù)),15/23順序(y,x---結(jié)果是先常數(shù)項,與polyfit相反)y為列向量;x為矩陣,第一列為全1列(即對應(yīng)于常數(shù)項),其余每一列對應(yīng)于一個變量(或一個含變量的項),即x要配成目標(biāo)函數(shù)的形式(常數(shù)項在最前)x中有多少列則結(jié)果的函數(shù)中就有多少項首先要確定要擬合的函數(shù)形式,然后確定待定的系,從常數(shù)項開始排列,須構(gòu)造x(每列對應(yīng)于函數(shù)中的一項,剔除待定系數(shù)),擬合就是確定待定系數(shù)的過程(當(dāng)然需先確定函數(shù)的型式)重點:regress(y,x)重點與難點是如何加工處理矩陣x。y是函數(shù)值,一定是只有一列。也即目標(biāo)函數(shù)的形式是由矩陣X來確定如s=a+b*x1+c*x2+d*x3+e*x1^2+f*x2*x3+g*x1^2,一定有一個常數(shù)項,且必須放在最前面(即x的第一列為全1列)X中的每一列對應(yīng)于目標(biāo)函數(shù)中的一項(目標(biāo)函數(shù)有多少項則x中就有多少列)X=[ones,x1,x2,x3,x1.^2,x2.*x3,x1.?2](剔除待定系數(shù)的形式)regress:y/x順序,矩陣X需要加工處理nlinfit:x/y順序,X/Y就是原始的數(shù)據(jù),不要做任何的加工。(即regress靠矩陣X來確定目標(biāo)函數(shù)的類型形式(所以X很復(fù)雜,要作很多處理)而nlinfit是靠程序來確定目標(biāo)函數(shù)的類型形式(所以X就是原始數(shù)據(jù),不要做任何處理)例1測16名成年女子的身高與腿長所得數(shù)據(jù)如下:身高143145146147150153154155156158159160162164腿長8885889192939395969897969899100102配成y=a+b*x形式>>x=[143145146147149150153154155156157158159160162164]';>>y=[8885889192939395969897969899100102]';>>plot(x,y,'r+')>>z=x;>>x=[ones(16,1),x];----常數(shù)項>>[b,bint,r,rint,stats]=regress(y,x);---處結(jié)果與polyfit(x,y,1)相同>>b,bint,stats得結(jié)果:b=bint=-16.0730-33.70711.5612------每一行為一個區(qū)間0.71940.60470.8340stats=0.9282180.95310.0000即;的置信區(qū)間為[-33.7017,1.5612],的置信區(qū)間為[0.6047,0.834];r2=0.9282,F=180.9531,p=0.0。p<0.05,可知回歸模型y=-16.073+0.7194x成立.>>[b,bint,r,rint,stats]=regress(Y,X,0.05);-----結(jié)果相同>>[b,bint,r,rint,stats]=regress(Y,X,0.03);>>polyfit(x,y,1)-----當(dāng)為一元時(也只有一組數(shù)),則結(jié)果與regress是相同的,只是命令中x,y要交換順序,結(jié)果的系數(shù)排列順序完全相反,x中不需要全1列。ans=0.7194-16.0730--此題也可用polyfit求解,殺雞用牛刀,脖子被切斷。3、殘差分析,作殘差圖:>>rcoplot(r,rint)從殘差圖可以看出,除第二個數(shù)據(jù)外,其余數(shù)據(jù)的殘差離零點均較近,且殘差的置信區(qū)間均包含零點,這說明回歸模型y=-16.073+0.7194x能較好的符合原始數(shù)據(jù),而第二個數(shù)據(jù)可視為異常點(而剔除)4、預(yù)測及作圖:>>plot(x,y,'r+')>>holdon>>a=140:165;>>b=b(1)+b(2)*a;>>plot(a,b,'g')例2觀測物體降落的距離s與時間t的關(guān)系,得到數(shù)據(jù)如下表,求s關(guān)于t的回歸方程t(s)1/302/303/304/305/306/307/30s(cm)11.8615.6720.6026.6933.7141.9351.13t(s)8/309/3010/3011/3012/3013/3014/30s(cm)61.4972.9085.4499.08113.77129.54146.48法一:直接作二次多項式回歸t=1/30:1/30:14/30;s=[11.8615.6720.6026.6933.7141.9351.1361.4972.9085.4499.08113.77129.54146.48];>>[p,S]=polyfit(t,s,2)p=489.294665.88969.1329得回歸模型為:方法二----化為多元線性回歸:t=1/30:1/30:14/30;s=[11.8615.6720.6026.6933.7141.9351.1361.4972.9085.4499.08113.77129.54146.48];>>T=[ones(14,1),t',(t.^2)']%???是否可行???等驗證...----因為有三個待定系數(shù),所以有三列,始于常數(shù)項>>[b,bint,r,rint,stats]=regress(s',T);>>b,statsb=9.132965.8896489.2946stats=1.0e+007*0.00001.037800.0000得回歸模型為:%結(jié)果與方法1相同>>T=[ones(14,1),t,(t.^2)']%???是否可行???等驗證...polyfit------一元多次regress----多元一次---其實通過技巧也可以多元多次regress最通用的,萬能的,表面上是多元一次,其實可以變?yōu)槎嘣啻吻胰我夂瘮?shù),如x有n列(不含全1列),則表達(dá)式中就有n+1列(第一個為常數(shù)項,其他每項與x的列序相對應(yīng))??????此處的說法需進一步驗證證……………例3設(shè)某商品的需求量與消費者的平均收入、商品價格的統(tǒng)計數(shù)據(jù)如下,建立回歸模型,預(yù)測平均收入為1000、價格為6時的商品需求量.需求量10075807050659010011060收入10006001200500300400130011001300300價格5766875439選擇純二次模型,即----用戶可以任意設(shè)計函數(shù)>>x1=[10006001200500300400130011001300300];>>x2=[5766875439];>>y=[10075807050659010011060]';X=[ones(10,1)x1'x2'(x1.^2)'(x2.^2)'];%注意技巧性??????[b,bint,r,rint,stats]=regress(y,X);%這是萬能方法??????需進一步驗證>>b,statsb=110.53130.1464-26.5709-0.00011.8475stats=0.970240.66560.000520.5771故回歸模型為:剩余標(biāo)準(zhǔn)差為4.5362,說明此回歸模型的顯著性較好.--------(此題還可以用rstool(X,Y)命令求解,詳見回歸問題詳解)>>X=[ones(10,1)x1'x2'(x1.^2)'(x2.^2)',sin(x1.*x2)',(x1.*exp(x2))'];>>[b,bint,r,rint,stats]=regress(y,X);>>b,stats(個人2011年認(rèn)為,regress只能用于函數(shù)中的每一項只能有一個待定系數(shù)的情況,不能用于aebx等的情況)regress(y,x)----re是y/x逆置的---y是列向量---須確定目標(biāo)函數(shù)的形式---x須構(gòu)造(通過構(gòu)造來反映目標(biāo)函數(shù))---x中的每一列與目標(biāo)函數(shù)的一項對應(yīng)(剔除待定系數(shù))----首項為常數(shù)項(x的第一列為全1)----有函數(shù)有n項(待定系數(shù)),則x就有n列----regress只能解決每項只有一個待定系數(shù)的情況且必須有常數(shù)項的情況(且每項只有一個待定系數(shù),即項數(shù)與待定系數(shù)數(shù)目相同)***其重(難、關(guān)鍵)點:列向量、構(gòu)造矩陣(X):目標(biāo)函數(shù)中的每項與X中的一列對應(yīng)。(由X來確定目標(biāo)函數(shù)的類型/形式)三、非線性回歸(擬合)使用格式:beta=nlinfit(x,y,‘程序名’,beta0)[beta,r,J]=nlinfit(X,y,fun,beta0)

X給定的自變量數(shù)據(jù),Y給定的因變量數(shù)據(jù),fun要擬合的函數(shù)模型(句柄函數(shù)或者內(nèi)聯(lián)函數(shù)形式),beta0函數(shù)模型中待定系數(shù)估計初值(即程序的初始實參)beta返回擬合后的待定系數(shù)其中beta為估計出的回歸系數(shù);r為殘差;J為Jacobian矩陣輸入數(shù)據(jù)x、y分別為n*m矩陣和n維列向量,對一元非線性回歸,x為n維列向量?!痬odel’為是事先用m-文件定義的非線性函數(shù);beta0為回歸系數(shù)的初值可以擬合成任意函數(shù)。最通用的,萬能的命令x,y順序,x不需要任何加工,直接用原始數(shù)據(jù)。(也不需要全1列)---所編的程序一定是兩個形參(待定系數(shù)/向量,自變量/矩陣:每一列為一個自變量)結(jié)果要看殘差的大小和是否有警告信息,如有警告則換一個b0初始向量再重新計算。本程序中也可能要用.*./.^如結(jié)果中有警告信息,則必須多次換初值來試算難點是編程序與初值nlinfit多元任意函數(shù),(自己任意設(shè)計函數(shù),再求待定系數(shù))順序([b,r,j]=nlinfit(x,y,’…’,b0)y為列向量;x為矩陣,無需加全1列,x,y就是原始的數(shù)據(jù)點,(x/y正順序,所以x不要加全1列)需預(yù)先編程(兩個參數(shù),系數(shù)向量,各變量的矩陣/每列為一個變量)存在的問題:不同的beta0,則會產(chǎn)生不同的結(jié)果,如何給待定系數(shù)的初值以及如何分析結(jié)果的好壞,如出現(xiàn)警告信息,則換一個待定系數(shù)試一試。因為擬合本來就是近似的,可能有多個結(jié)果。1:重點(難點)是預(yù)先編程序(即確定目標(biāo)函數(shù)的形式,而regress的目標(biāo)函數(shù)由x矩陣來確定,其重難點為構(gòu)造矩陣a)2:x/y順序—列向量----x/y是原始數(shù)據(jù),不要做任何修改3:編程:一定兩個形參(beta,x)a=beta(1);b=beta(2);c=beta(3);…x1=x(:,1);x2=x(:,2);x3=x(:,3);即每一列為一個自變量4:regress/nlinfit都是列向量5:regress:有n項(n個待定系數(shù)),x就有n列;nlinfit:有m個變量則x就有m列例1已知數(shù)據(jù):x1=[0.5,0.4,0.3,0.2,0.1];

x2=[0.3,0.5,0.2,0.4,0.6];x3=[1.8,1.4,1.0,1.4,1.8];y=[0.785,0.703,0.583,0.571,0.126]’;且y與x1,x2,x3關(guān)系為多元非線性關(guān)系(只與x2,x3相關(guān))為:

y=a+b*x2+c*x3+d*(x2.^2)+e*(x3.^2)—此函數(shù)是由用戶根據(jù)圖形的形狀等所配的曲線,即自己選定函數(shù)類型求非線性回歸系數(shù)a,b,c,d,e。(1)對回歸模型建立M文件model.m如下:functionyy=myfun(beta,x)%一定是兩個參數(shù):系數(shù)和自變量---一個向量/一個矩陣a=beta(1)b=beta(2)c=beta(3)x1=x(:,1);%系數(shù)是數(shù)組,b(1),b(2),…b(n)依次代表系數(shù)1,系數(shù)2,……系數(shù)nx2=x(:,2);%自變量x是一個矩陣,它的每一列分別代表一個變量,有n列就可以最多nx3=x(:,3);yy=beta(1)+beta(2)*x2+beta(3)*x3+beta(4)*(x2.^2)+beta(5)*(x3.^2);(b(i)與待定系數(shù)的順序關(guān)系可以任意排列,并不是一定常數(shù)項在最前,只是結(jié)果與自己指定的相對應(yīng))(x一定是一列對應(yīng)一個變量,不能x1=x(1),x2=x(2),x3=x(3)……)(2)主程序如下:x=[0.5,0.4,0.3,0.2,0.1;0.3,0.5,0.2,0.4,0.6;1.8,1.4,1.0,1.4,1.8]';-----每一列為一個變量y=[0.785,0.703,0.583,0.571,0.126]';beta0=[1,1,1,1,1,1]';%有多少個待定系數(shù),就給多少個初始值。[beta,r,j]=nlinfit(x,y,@myfun,beta0)beta=-0.44205.51110.3837-8.1734-0.1340此題也可用regress來求解,但結(jié)果是不一樣的>>x1=[0.5,0.4,0.3,0.2,0.1];>>x2=[0.3,0.5,0.2,0.4,0.6];>>x3=[1.8,1.4,1.0,1.4,1.8];>>y=[0.785,0.703,0.583,0.571,0.126]';>>n=length(x1);>>x=[ones(n,1),x2',x3',(x2.^2)',(x3.^2)'];>>[b,bint,r,rint,stats]=regress(y,x);>>b,statsb=-3.3844-1.84506.51370-2.1773stats=0.78591.22320.56740.05572011年題目改為:y=a+b*x1+c*x2+d*(x3.^2)+e*(x1.^2)+f*sin(x2)求非線性回歸系數(shù)a,b,c,d,e,ffunctionf=fxxnh(beta,x)%所編的程序一定是兩個形參,第一個為待定系數(shù)向量,第二個為自變量矩陣a=beta(1);b=beta(2);c=beta(3);d=beta(4);e=beta(5);f=beta(6);%系數(shù)向量中的一個元素代表一個待定系數(shù)x1=x(:,1);%自變量矩陣每一列代表一個自變量x2=x(:,2);x3=x(:,3);f=a+b.*x1+c.*x2+d.*(x3.^2)+e.*(x1.^2)+f.*sin(x2);但計算出現(xiàn)了問題例2混凝土的抗壓強度隨養(yǎng)護時間的延長而增加,現(xiàn)將一批混凝土作成12個試塊,記錄了養(yǎng)護日期(日)及抗壓強度y(kg/cm2)的數(shù)據(jù):養(yǎng)護時間:x=[234579121417212856]抗壓強度:y=[35+r42+r47+r53+r59+r65+r68+r73+r76+r82+r86+r99+r]建立非線性回歸模型,對得到的模型和系數(shù)進行檢驗。注明:此題中的+r代表加上一個[-0.5,0.5]之間的隨機數(shù)模型為:y=a+k1*exp(m*x)+k2*exp(-m*x);------有四個待定系數(shù)Matlab程序:

x=[234579121417212856];

r=rand(1,12)-0.5;

y1=[354247535965687376828699];

y=y1+r;

myfunc=inline('beta(1)+beta(2)*exp(beta(4)*x)+beta(3)*exp(-beta(4)*x)','beta','x');-----beta=nlinfit(x,y,myfunc,[0.50.50.50.5]);---初值為0.2也可以,如為1則不行,則試著換系數(shù)初值----此處為一元,x’,y’行/列向量都可以a=beta(1),k1=beta(2),k2=beta(3),m=beta(4)%testthemodelxx=min(x):max(x);-----2:56

yy=a+k1*exp(m*xx)+k2*exp(-m*xx);

plot(x,y,'o',xx,yy,'r')

結(jié)果:

a=87.5244

k1=0.0269

k2=-63.4591

m=0.1083

圖形:此題不能用regress求解,因為有些式子中含有兩個待定系數(shù)例3出鋼時所用的盛鋼水的鋼包,由于鋼水對耐火材料的侵蝕,容積不斷增大.我們希望知道使用次數(shù)與增大的容積之間的關(guān)系.對一鋼包作試驗,測得的數(shù)據(jù)列于下表:使用次數(shù)增大容積使用次數(shù)增大容積234567896.428.209.589.509.7010.009.939.991011121314151610.4910.5910.6010.8010.6010.9010.76對將要擬合的非線性模型y=aeb/x,(如再加y=c*sin(x)+aeb/x)建立m-文件volum.m如下:functionyhat=volum(beta,x)yhat=beta(1)*exp(beta(2)./x);或functionf=zhang1(beta,x)a=beta(1);b=beta(2);f=a*exp(b./x);---2、輸入數(shù)據(jù):>>x=2:16;>>y=[6.428.209.589.59.7109.939.9910.4910.5910.6010.8010.6010.9010.76];>>beta0=[82]';----初值[1,1]也可以3、求回歸系數(shù):>>[beta,r,J]=nlinfit(x',y','volum',beta0);%beta0初值為列/行向量都可以,還是為列吧。>>betabeta=11.6037-1.0641即得回歸模型為:4、預(yù)測及作圖:>>[YY,delta]=nlpredci('volum',x',beta,r,J)>>plot(x,y,'k+',x,YY,'r')或>>plot(x,y,'ro')>>holdon>>xx=2:0.05:16;>>yy=beta(1)*exp(beta(2)./xx);>>plot(xx,yy,'g')又或>>plot(x,y,'ro')>>holdon>>xx=2:0.05:16;>>yy=volum(beta,xx);--------通過調(diào)用用戶自編的函數(shù)>>plot(xx,yy,'g')>>[beta,r,J]=nlinfit(x',y','volum',[1,1]);%下面換了多個初值,結(jié)果都是一樣的。>>betabeta=11.6037-1.0641>>[beta,r,J]=nlinfit(x',y','volum',[1,5]);>>betabeta=11.6037-1.064>>[beta,r,J]=nlinfit(x',y','volum',[10,5]);beta=11.6037-1.0641>>[beta,r,J]=nlinfit(x',y','volum',[10,50]);beta=11.6037-1.0641以下用來lsqcurvefit求解,結(jié)果是一樣的。>>[beta,a,b,exitflag]=lsqcurvefit('volum',[8,2],x',y')Optimizationterminated:relativefunctionvaluechangingbylessthanOPTIONS.TolFun.beta=11.6037-1.0641exitflag=3>>[beta,a,b,exitflag]=lsqcurvefit('volum',[1,1],x',y')%換不同的初值,結(jié)果是一樣的。beta=11.6037-1.0641exitflag=3>>[beta,a,b,exitflag]=lsqcurvefit('volum',[10,1],x',y')beta=11.6037-1.0641exitflag=3>>[beta,a,b,exitflag]=lsqcurvefit('volum',[10,5],x',y')beta=11.6037-1.0641exitflag=3>>[beta,a,b,exitflag]=lsqcurvefit('volum',[10,50],x',y')beta=11.6037-1.0641exitflag=3例4財政收入預(yù)測問題:財政收入與國民收入、工業(yè)總產(chǎn)值、農(nóng)業(yè)總產(chǎn)值、總?cè)丝?、就業(yè)人口、固定資產(chǎn)投資等因素有關(guān)。下表列出了1952-1981年的原始數(shù)據(jù),試構(gòu)造預(yù)測模型。財政收入預(yù)測問題:財政收入與國民收入、工業(yè)總產(chǎn)值、農(nóng)業(yè)總產(chǎn)值、總?cè)丝凇⒕蜆I(yè)人口、固定資產(chǎn)投資等因素有關(guān)。下表列出了1952-1981年的原始數(shù)據(jù),試構(gòu)造預(yù)測模型。年份國民收入(億元)工業(yè)總產(chǎn)值(億元)農(nóng)業(yè)總產(chǎn)值(億元)總?cè)丝冢ㄈf人)就業(yè)人口(萬人)固定資產(chǎn)投資(億元)財政收入(億元)195259834946157482207294418419535864554755879621364892161954707520491602662183297248195573755852961465223289825419568257155566282823018150268195783779857564653237111392861958102812355986599426600256357195911141681509672072617333844419601079187044466207258803805061961757115643465859255901382711962677964461672952511066230196377910465146917226640852661964943125058470499277361293231965115215816327253828670175393196613221911687745422980521246619671249164769776368308141563521968118715656807853431915127303196913722101688806713322520744719701638274776782992344323125641971178031567908522935620355638197218333365789871773585435465819731978368485589211366523746911974199336968919085937369393655197521214254932924213816846269219762052430995593717388344436571977218949259719497439377454723197824755590105896259398565509221979270260651150975424058156489019802791659211949870541896568826198129276862127310007273280496810解設(shè)國民收入、工業(yè)總產(chǎn)值、農(nóng)業(yè)總產(chǎn)值、總?cè)丝?、就業(yè)人口、固定資產(chǎn)投資分別為x1、x2、x3、x4、x5、x6,財政收入為y,設(shè)變量之間的關(guān)系為:y=ax1+bx2+cx3+dx4+ex5+fx6使用非線性回歸方法求解。1.

對回歸模型建立M文件model.m如下:functionyy=model(beta0,X)%一定是兩個參數(shù),第一個為系數(shù)數(shù)組,b(1),b(2),…b(n)%分別代表每個系數(shù),而第二個參數(shù)代表所有的自變量,%是一個矩陣,它的每一列分別代表一個自變量。a=beta0(1);b=beta0(2);%每個元素c=beta0(3);d=beta0(4);e=beta0(5);f=beta0(6);x1=X(:,1);%每一列x2=X(:,2);x3=X(:,3);x4=X(:,4);x5=X(:,5);x6=X(:,6);yy=a*x1+b*x2+c*x3+d*x4+e*x5+f*x6;2.主程序liti6.m如下:X=[598.00,349.00,461.00,57482.00,20729.00,44.00;586,455,475,58796,21364,89;707,520,491,60266,21832,97;737,558,529,61465,22328,98;825,715,556,62828,23018,150;837,798,575,64653,23711,139;1028,1235,598,65994,26600,256;1114,1681,509,67207,26173,338;1079,1870,444,66207,25880,380;757,1156,434,65859,25590,138;677,964,461,67295,25110,66;779,1046,514,69172,26640,85;943,1250,584,70499,27736,129;1152,1581,632,72538,28670,175;1322,1911,687,74542,29805,212;1249,1647,697,76368,30814,156;1187,1565,680,78534,31915,127;1372,2101,688,80671,33225,207;1638,2747,767,82992,34432,312;1780,3156,790,85229,35620,355;1833,3365,789,87177,35854,354;1978,3684,855,89211,36652,374;1993,3696,891,90859,37369,393;2121,4254,932,92421,38168,462;2052,4309,955,93717,38834,443;2189,4925,971,94974,39377,454;2475,5590,1058,96259,39856,550;2702,6065,1150,97542,40581,564;2791,6592,1194,98705,41896,568;2927,6862,1273,100072,73280,496];y=[184.00216.00248.00254.00268.00286.00357.00444.00506.00...271.00230.00266.00323.00393.00466.00352.00303.00447.00...564.00638.00658.00691.00655.00692.00657.00723.00922.00...890.00826.00810.0]';beta0=[0.50-0.03-0.600.01-0.020.35];betafit=nlinfit(X,y,'model',beta0)結(jié)果為betafit=0.5243-0.0294-0.63040.0112-0.02300.3658(結(jié)果也可能是:0.3459-0.0180-0.37000.0030-0.00200.4728)即y=0.5243x1-0.0294x2-0.6304x3+0.0112x4-0.0230x5+0.3658x6此題也可以用regress來求解(我自己做的,不一定對???)----結(jié)果有些不同,含有一個常數(shù)>>clear>>x=xlsread('cz.xls');%已經(jīng)把所有的有效數(shù)據(jù)拷入到cd.xls文件中去了。>>y=x(:,7);>>x(:,7)=[];>>z=ones(30,1);>>x=[z,x];>>[b,bint,r,rint,states]=regress(y,x);>>b,statesb=159.14400.4585-0.0112-0.51250.0008-0.00280.3165stats=1.0e+003*0.00100.228301.0488X,y的原始數(shù)據(jù):helpnlinfit/helpnlinfitExamples:Use@tospecifyMODELFUN:loadreaction;beta=nlinfit(reactants,rate,@mymodel,beta);whereMYMODELisaMATLABfunctionsuchas:functionyhat=mymodel(beta,x)yhat=(beta(1)*x(:,2)-x(:,3)/beta(5))./...(1+beta(2)*x(:,1)+beta(3)*x(:,2)+beta(4)*x(:,3));ExamplesFUNcanbespecifiedusing@:nlintool(x,y,@myfun,b0)whereMYFUNisaMATLABfunctionsuchas:functionyhat=myfun(beta,x)b1=beta(1);b2=beta(2);yhat=1./(1+exp(b1+b2*x));fun=inline('1./(1+exp(b(1)+b(2)*x))','b','x')nlintool(x,y,fun,b0)每一列為一個變量,所以有n個自變量,就有n列X1,x2,…xn,y都為列向量四、非線性回歸或曲線回歸問題配曲線的一般方法是:(一)先對兩個變量x和y作n次試驗觀察得畫出散點圖,散點圖(二)根據(jù)散點圖確定須配曲線的類型.通常選擇的六類曲線如下:(1)雙曲線(2)冪函數(shù)曲線y=a,其中x>0,a>0(3)指數(shù)曲線y=a其中參數(shù)a>0.(4)倒指數(shù)曲線y=a其中a>0,(5)對數(shù)曲線y=a+blogx,x>0(6)S型曲線(三)然后由n對試驗數(shù)據(jù)確定每一類曲線的未知參數(shù)a和b.解例2.由散點圖我們選配倒指數(shù)曲線y=a根據(jù)線性化方法,算得由此最后得作業(yè)1、考察溫度x對產(chǎn)量y的影響,測得下列10組數(shù)據(jù):溫度(℃)20253035404550556065產(chǎn)量(kg)17.117.918.719.621.222.524.3求y關(guān)于x的線性回歸方程,檢驗回歸效果是否顯著,并預(yù)測x=42℃時產(chǎn)量的估值及預(yù)測區(qū)間(置信度95%).2、某零件上有一段曲線,為了在程序控制機床上加工這一零件,需要求這段曲線的解析表達(dá)式,在曲線橫坐標(biāo)xi處測得縱坐標(biāo)yi共11對數(shù)據(jù)如下:xi02468101214161820yi0.62.04.47.511.817.123.331.239.649.761.7求這段曲線的縱坐標(biāo)y關(guān)于橫坐標(biāo)x的二次多項式回歸方程.3:在研究化學(xué)動力學(xué)反應(yīng)過程中,建立了一個反應(yīng)速度和反應(yīng)物含量的數(shù)學(xué)模型,形式為其中是未知參數(shù),是三種反應(yīng)物(氫,n戊烷,異構(gòu)戊烷)的含量,y是反應(yīng)速度.今測得一組數(shù)據(jù)如表4,試由此確定參數(shù),并給出置信區(qū)間.的參考值為(1,0.05,0.02,0.1,2)序號反應(yīng)速度y氫x1n戊烷x2異構(gòu)戊烷x318.554703001023.79285801034.8247030012040.024708012052.754708010614.391001901072.54100806584.3547019065913.0010030054108.50100300120110.05100801201211.3228530010133.132851901204、混凝土的抗壓強度隨養(yǎng)護時間的延長而增加,現(xiàn)將一批混凝土作成12個試塊,記錄了養(yǎng)護日期x(日)及抗壓強度y(kg/cm2)的數(shù)據(jù):養(yǎng)護時間x234579121417212856抗壓強度y354247535965687376828699試求型回歸方程.關(guān)于利用Matab進行指定非線性擬合的問題(轉(zhuǎn)載)

(2010-05-2716:46:07)轉(zhuǎn)載▼標(biāo)簽:

雜談分類:

學(xué)術(shù)※1。優(yōu)化工具箱的利用

函數(shù)

描述

LSQLIN

有約束線性最小二乘優(yōu)化

LSQNONNEG

非負(fù)約束線性最小二乘優(yōu)化問題

當(dāng)有約束問題存在的時候,應(yīng)該采用上面的方法代替Polyfit與反斜線(\)。具體例子請參閱優(yōu)化工具箱文檔中的相應(yīng)利用這兩個函數(shù)的例子。d.

非線性曲線擬合

利用MATLAB的內(nèi)建函數(shù)

函數(shù)名

描述

FMINBND

只解決單變量固定區(qū)域的最小值問題

FMINSEARCH多變量無約束非線性最小化問題(Nelder-Mead方法)。

下面給出一個小例子展示一下如何利用FMINSEARCH

1.首先生成數(shù)據(jù)

>>t=0:.1:10;

>>t=t(:);

>>Data=40*exp(-.5*t)+rand(size(t));

%將數(shù)據(jù)加上隨機噪聲

2.寫一個m文件,以曲線參數(shù)作為輸入,以擬合誤差作為輸出

functionsse=myfit(params,Input,Actural_Output)

A=params(1);

lamda=params(2);

Fitted_Curve=A.*exp(-lamda*Input);

Error_Vector=Fitted_Curve-Actural_Output;

%當(dāng)曲線擬合的時候,一個典型的質(zhì)量評價標(biāo)準(zhǔn)就是誤差平方和

sse=sum(Error_Vector.^2);

%當(dāng)然,也可以將sse寫作:sse=Error_Vector(:)*Error_Vector(:);3.調(diào)用FMINSEARCH

>>Strarting=rand(1,2);

>>options=optimset('Display','iter');

>>Estimates=fiminsearch(@myfit,Strarting,options,t,Data);

>>plot(t,Data,'*');

>>holdon

>>plot(t,Estimates(1)*exp(-Estimates(2)*t),'r');

Estimates將是一個包含了對原數(shù)據(jù)集進行估計的參數(shù)值的向量。

附圖見后面:

FMINSEARCH通常能夠用來解決不連續(xù)情況,特別是如果他們不出現(xiàn)在解的附近的時候。它得到的通常也是局部解。FMINSEARCH只能夠最小化實數(shù)值(也就是說,解的域必須只能包括實數(shù),函數(shù)的輸出只能夠為實數(shù)值)。當(dāng)感興趣的是復(fù)數(shù)變量的域的時候,他們必須被分割為實部與虛部?!?.MATLAB的FIGURE窗口:最基本的擬合界面與數(shù)據(jù)統(tǒng)計工具

MATLAB通過基本的擬合界面也支持基本曲線擬合。利用這個界面,你可以快速地在簡單易用的環(huán)境中實現(xiàn)許多基本的曲線擬合。這個界面可以實現(xiàn)以下功能:

a.通過比樣條插值(splineinterpolant)、hermite插值、或者是高達(dá)10階的多項式插值實現(xiàn)數(shù)據(jù)的擬合;

b.對給定數(shù)據(jù)同時實現(xiàn)多樣插值的繪制;

c.繪制殘差圖;

d.檢查擬合結(jié)果的殘差的數(shù)值;

e.通過內(nèi)插值或者外推插值評價一個擬合結(jié)果;

f.對擬合結(jié)果和殘差的模進行圖形繪制;

g.將擬合結(jié)果保存入MATLAB工作空間。

開發(fā)你的擬合應(yīng)用的時候,你可以通過基本擬合(BasicFitting)界面,也可以通過命令行函數(shù),也可以同時作用。你可以通過基本擬合界面只能夠?qū)崿F(xiàn)2-D數(shù)據(jù)的擬合。然而,如果你用subplot繪制多個數(shù)據(jù)集,只要有至少一個數(shù)據(jù)集是2D的,那么就可以用基本擬合界面。

可以通過如下步驟激活基本擬合界面:

1.繪制數(shù)據(jù);

2.從figure窗口的Tools菜單條下面選擇BasicFitting菜單項;

有關(guān)BasicFitting界面的更多信息,請查閱MATLAB幫助文檔的相應(yīng)部分。

注意:對于HP,IBM以及SGI平臺,MATLAB6.0(R12.0)以及MATLAB6.1(R12.1)的基本擬合界面不受支持。

數(shù)據(jù)統(tǒng)計界面可以用來對圖形中的每個數(shù)據(jù)集進行統(tǒng)計量的計算??梢酝ㄟ^如下步驟將數(shù)據(jù)統(tǒng)計界面激活:

1.制數(shù)據(jù);

2.從figure窗口的Tools菜單條下面選擇DataStatistics菜單項;

一。優(yōu)化工具箱函數(shù)

LSQNONLIN解決非線性最小二乘法問題,包括非線性數(shù)據(jù)擬合問題

LSQCURVEFIT

解決非線性數(shù)據(jù)擬合問題

下面給出利用這兩個函數(shù)的例子:

LSQNONLIN:利用這個函數(shù)最小化連續(xù)函數(shù)只能夠找到句柄解。下面的例子說明利用LSQNONLIN函數(shù)用下面的函數(shù)進行擬合:

f=A+Bexp(C*x)+D*exp(E*x)

對數(shù)據(jù)集x與y進行擬合,其中y是在給定x的情況下的期望輸出(可以是方程給出數(shù)組,也可以是單獨數(shù)據(jù)組成的數(shù)組)。為了解決這個問題,先建立下面的名為fit_simp.m的函數(shù),它利用數(shù)據(jù)x與y,將他們作為優(yōu)化輸入?yún)?shù)傳遞給LSQNONLIN。利用給定的數(shù)據(jù)x計算f的值,再與原始數(shù)據(jù)y進行比較。經(jīng)驗值與實際計算出的值之間的差異作為輸出值返回。LSQNOLIN函數(shù)就是最小化這些差的平方和。

function

diff=fit_simp(x,X,Y)

%此函數(shù)被LSQNONLIN調(diào)用

%x是包含等式系數(shù)的向量

%X與Y是作為操作數(shù)傳遞給lsnonlin

A=x(1);

B=x(2);

C=x(3);

D=x(4);

E=x(5);

diff=A+B.*exp(C.*X)+D.*exp(E.*X)-Y;下面的腳本是利用上面定義的fit_simp.m函數(shù)的一個例子:

%定義你打算擬合的數(shù)據(jù)集合

>>X=0:.01:.5;

>>Y=2.0.*exp(5.0.*X)+3.0.*exp(2.5.*X)+1.5.*rand(size(X));

%初始化方程系數(shù)

>>X0=[11111]';

%設(shè)置用中等模式(memdium-scale)算法

>>options=optimset('Largescale','off');

%通過調(diào)用LSQNONLIN重現(xiàn)計算新的系數(shù)

>>x=lsqnonlin(@fit_simp,X0,[],[],options,X,Y);

%調(diào)用LSQNONLIN結(jié)果輸出表明擬合是成功的

Optimizationterminatedsuccessfully:

GradientinthesearchdirectionlessthantolFun

Gradientlessthan10*(tolFun+tolX)

%繪制原始數(shù)據(jù)與新的計算的數(shù)據(jù)

>>Y_new=x(1)+x(2).*exp(x(3).*X)+x(4).*exp(x(5).*X);

>>plot(X,Y,'+r',X,Y_new,'b');※注意:LSQNONLIN只可以處理實數(shù)變量。在處理包括復(fù)數(shù)變量的實例的擬合的時候,數(shù)據(jù)集應(yīng)該被切分成實數(shù)與虛數(shù)部分。下面給出一個例子演示如何對復(fù)數(shù)參數(shù)進行最小二乘擬合。

為了擬合復(fù)數(shù)變量,你需要將復(fù)數(shù)分解為實數(shù)部分與虛數(shù)部分,然后把他們傳遞到函數(shù)中去,這個函數(shù)被LEASTSQ作為單個輸入調(diào)用。首先,將復(fù)數(shù)分解為實部與虛部兩個向量。其次,將這兩個向量理解成諸如第一部分是實部、第二部分是虛部。在MATLAB函數(shù)中,重新裝配復(fù)數(shù)數(shù)據(jù),并用你想擬合的復(fù)數(shù)方程計算。將輸出向量分解實部與虛部,將這兩部分連接為一個單一的輸出向量傳遞回LEASTSQ。下面,給出一個例子演示如何根據(jù)兩個復(fù)數(shù)指數(shù)擬合實數(shù)X與Y。

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會有圖紙預(yù)覽,若沒有圖紙預(yù)覽就沒有圖紙。
  • 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
  • 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護處理,對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對任何下載內(nèi)容負(fù)責(zé)。
  • 6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時也不承擔(dān)用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

評論

0/150

提交評論