Matlab解決線性回歸問題論文

Matlab解決線性回歸問題2010級(jí)電子信息工程班【摘要】MATLAB和Mathematica、Maple并稱為三大數(shù)學(xué)軟件。它在數(shù)學(xué)類科技應(yīng)用軟件中在數(shù)值計(jì)算方面首屈一指°MATLAB可以進(jìn)行矩陣運(yùn)算、繪制函數(shù)和數(shù)據(jù)、實(shí)現(xiàn)算法、創(chuàng)建用戶界面、連接其他編程語言的程序等，主要應(yīng)用于工程計(jì)算、控制設(shè)計(jì)、信號(hào)處理與通訊、圖像處理、信號(hào)檢測(cè)、金融建模設(shè)計(jì)與分析等領(lǐng)域?；貧w分析，是對(duì)現(xiàn)有數(shù)據(jù)進(jìn)行處理、從中發(fā)現(xiàn)有用信息的一種重要手段。而線性回歸，特別是一元線性回歸分析更是人們優(yōu)先考慮采用的方式?；诖?，本文就一元線性回歸的MATLAB實(shí)現(xiàn)作了一番探討，給出了多種實(shí)現(xiàn)方式，并通過一個(gè)實(shí)例加以具體展示，在數(shù)據(jù)處理時(shí)可根據(jù)自己的需要靈活地加以選用?！娟P(guān)鍵詞】MATLAB程序、一元線性回歸、多元線性回歸問題、線性回歸擬合問題、Subplot、三次樣條插值函數(shù)。一、 提出問題在回歸分析中，如果有兩個(gè)或兩個(gè)以上的自變量，就稱為多元回歸。事實(shí)上，一種現(xiàn)象常常是與多個(gè)因素相聯(lián)系的，由多個(gè)自變量的最優(yōu)組合共同來預(yù)測(cè)或估計(jì)因變量，比只用一個(gè)自變量進(jìn)行預(yù)測(cè)或估計(jì)更有效，更符合實(shí)際。二、 簡(jiǎn)述原理（一）一元線性回歸命令polyfit最小二乘多項(xiàng)式擬合[p，S]=polyfit（x，y，m）多項(xiàng)式y(tǒng)二a1xm+a2xm—1+…+amx+am+1其中x二（xl,x2,…，xm）xl???xm為（n*l）的矩陣；y為（n*l）的矩陣；p=（al,a2,…，am+1）是多項(xiàng)式y(tǒng)二a1xm+a2xm-1+???+amx+am+1的系數(shù)；S是一個(gè)矩陣，用來估計(jì)預(yù)測(cè)誤差.命令polyval多項(xiàng)式函數(shù)的預(yù)測(cè)值Y二polyval（p,x）求polyfit所得的回歸多項(xiàng)式在x處的預(yù)測(cè)值Y；p是polyfit函數(shù)的返回值；x和polyfit函數(shù)的x值相同。命令polyconf殘差個(gè)案次序圖[Y,DELTA]二polyconf（p,x,S,alpha）求polyfit所得的回歸多項(xiàng)式在x處的預(yù)測(cè)值Y及預(yù)測(cè)值的顯著性為1-alpha的置信區(qū)間DELTA；alpha缺省時(shí)為0.05。p是polyfit函數(shù)的返回值；x和polyfit函數(shù)的x值相同；S和polyfit函數(shù)的S值相同。命令polytool（x,y,m）—元多項(xiàng)式回歸命令命令regress多元線性回歸（可用于一元線性回歸）b=regress（Y,X）[b,bint,r,rint,stats]=regress（Y,X,alpha）b回歸系數(shù)bint回歸系數(shù)的區(qū)間估計(jì)r殘差rint殘差置信區(qū)間stats用于檢驗(yàn)回歸模型的統(tǒng)計(jì)量，有三個(gè)數(shù)值：相關(guān)系數(shù)R2、F值、與F對(duì)應(yīng)的概率p,相關(guān)系數(shù)R2越接近1,說明回歸方程越顯著；F>F1-a（k,n-k-1）時(shí)拒絕HO,F越大，說明回歸方程越顯著；與F對(duì)應(yīng)的概率p時(shí)拒絕HO,回歸模型成立。Y為n*1的矩陣；乂為（ones（n,1）,x1,…，xm）的矩陣；alpha顯著性水平（缺省時(shí)為0.05）。（二）多元線性回歸命令regress命令rstool多元二項(xiàng)式回歸命令：rstool（x,y,'model',alpha）x為n*m矩陣y為n維列向量model由下列4個(gè)模型中選擇1個(gè)（用字符串輸入,缺省時(shí)為線性模型）：linear（線性）：purequadratic（純二次）：interaction（交叉）：quadratic（完全二次）：alpha顯著性水平(缺省時(shí)為0.05)返回值beta系數(shù)返回值rmse剩余標(biāo)準(zhǔn)差返回值residuals殘差非線性回歸1．命令nlinfit[beta,R,J]=nlinfit(X,Y,''model',beta0)X為n*m矩陣Y為n維列向量model為自定義函數(shù)betaO為估計(jì)的模型系數(shù)beta為回歸系數(shù)R為殘差2．命令nlintoolnlintool(X,Y,'model',betaO,alpha)X為n*m矩陣Y為n維列向量model為自定義函數(shù)beta0為估計(jì)的模型系數(shù)alpha顯著性水平(缺省時(shí)為0.05)3．命令nlparcibetaci=nlparci(beta,R,J)beta為回歸系數(shù)R為殘差返回值為回歸系數(shù)beta的置信區(qū)間4．命令nlpredci[Y,DELTA]=nlpredci(‘model',X,beta,R,J)Y為預(yù)測(cè)值DELTA為預(yù)測(cè)值的顯著性為1-alpha的置信區(qū)間；alpha缺省時(shí)為0.05。X為n*m矩陣model為自定義函數(shù)beta為回歸系數(shù)R為殘差三、多元線性回歸問題在實(shí)際經(jīng)濟(jì)問題中，一個(gè)變量往往受到多個(gè)變量的影響。例如，家庭消費(fèi)支出，除了受家庭可支配收入的影響外，還受諸如家庭所有的財(cái)富、物價(jià)水平、金融機(jī)構(gòu)存款利息等多種因素的影響，表現(xiàn)在線性回歸模型中的解釋變量有多個(gè)。這樣的模型被稱為多元線性回歸模型。(multivariablelinearregressionmodel)1、多元線性回歸模型的一般形式多元線性回歸模型的一般形式為：

TOC\o"1-5"\h\zY=0+lXli+2X2i++kXki+i,i=l,2,,n (1)其中k為解釋變量的數(shù)目，j(jl,2,k)稱為回歸系(regressioncoefficient)。式也被稱為總體回歸函數(shù)的隨機(jī)表達(dá)式。它的非隨機(jī)表達(dá)式為：Y=0+1X1i+2X2i++kXki,i=1,2,n (2)j也被稱為偏回歸系數(shù)(partialregressioncoefficient)。2、多元線性回歸計(jì)算模型Y=0+lxl+2x2+kxk+N(0,2) (3)多元性回歸模型的參數(shù)估計(jì)，同一元線性回歸方程一樣，也是在要求誤差平方和(Ze)為最小的前提下，用最小二乘法或最大似然估計(jì)法求解參數(shù)。設(shè)(xll, x12 , x1 p, y1 ),,( xn1, xn2,xnp,yn)是一個(gè)樣本，用最大似然估計(jì)法估計(jì)參數(shù)：取b0,blbp，當(dāng)bObO,b1b1,bb達(dá)到最小。例1輸入以下程序:%數(shù)據(jù)y=[8SJ85.08,91,92.93,93,95.96,98,97.96,98,99,100,102]:3=[onesIlength(y)jl)js']%數(shù)據(jù)[b,bintjr：rint,stata]=reEre53(Y,X); %制圖輸出結(jié)果如圖1所示:

圖1圖1在Matlab命令窗口輸入z=b(l)+b(2)*x%運(yùn)算%輸出運(yùn)算plot(x,Y,‘k+，,x,z,，r，)%輸出運(yùn)算得到預(yù)測(cè)比較圖如圖2、3所示:

例2：（多元線性回歸）水泥凝固時(shí)放出的熱量y與水泥中的四種化學(xué)成分x1,x2,x3,x4有關(guān),今測(cè)得一組數(shù)據(jù)如下，試確定多元線性模型

序號(hào)盧和；冠8'、fip初；11^'P；1V'愛：也n2闘2沖亦5V71Px3^6^15P8r鬧存階込47V：33p6^7?,5^74.3^104.3^87.6^953109.2^1睦.斑序號(hào)戶扣伽他12*^T口公「21^侗11+J伽3佃珈4衍斗W盹x3^耳1腳4+1湖■:22卩26^珈12^72.5P93.143115.^83.8^109.4+1□分折：詛U413,1坯4,11J0J4J竝=[2伉29,托」1.52,557131,54,47,40^,15盯汁^=[615^8 1了衛(wèi),1S,4罟兌町屮材=[150,52^,4733^2,6,44^2益』4,12,12]屮尸卩8.5,743,104.3,87.6,95.9,109：2,102.77^->,93.1,115.9,8^.8,113.3,109.4]；^由武?應(yīng)可得X=[@TRT,J?T,0T護(hù)T](旳為單位列向星1屮Y=jT..Matlab程序?yàn)?:，:*■'d輸入勿]區(qū)命令：3小卩m?；廠門￡ %數(shù)據(jù)必呼護(hù)㈠“沐壬門“沖尹汽％數(shù)據(jù)J、：、?，?:.-■■-■二：■..■；'■ %數(shù)據(jù)■--■■：.二-'二■■-、…:-■=■- %數(shù)據(jù)y=[78.5.74.3,104.3.87.6.95.9J0^2.10^.7.72.5.93.U1S.9.&3.&.113.3.1014]>%數(shù)據(jù)X=[ones(length(y)J):x1,.x2\xJ\x4l：%把行向蛍轉(zhuǎn)鐵為列向魚屮丫二丫：％把行向HfeO列向量屮[tirtiintnr:rintzstatsj=regress“肉屮 %運(yùn)算公式b:bint;stats+J在Matlab圖示所示：心'. .'. ' ％數(shù)據(jù)X率[施29.55..31>52,55,71.丸51,47,迪66,68];灶&盡&血&9,17,22,也嗨網(wǎng);x4=[60,52,絢卸33,22禺咼22,繩軸12,121.注厲盡他3謝餌肌缶恫嗾嘰型疏人隆氐0$】昇遽瘀閱禺1注迢呃4];X=[onesQwifthfyh1）注1’注2*,訝注4’];%數(shù)據(jù)%數(shù)據(jù)%數(shù)據(jù)%運(yùn)行方程輸出結(jié)果如圖4所示:=62.40S4:p 的置信區(qū)間(t99.1786,2235893)b-.=1.5511^￡的置信區(qū)間(10.1663,3.26E5)/因此我們可得磯=05102」bnz的置信區(qū)間(11.1589,2.1792)^^=0.1019>陰的置信區(qū)間<n.63S5:1.S423)鋭=口441彳獲的置信區(qū)間017791：1.4910>口=0朋24』=111.4792^=0.0000.^p<W回歸模型」162.4054+1.5511x1+0.5102x2+0.1019^3-0.1441^立一結(jié)果如圖5所示:

圖5在matlab命令框中輸入衛(wèi)二feCHb僅門1+b(3)*x2+b(4}^3+t)譯片4屮plotfX.Y.-k+rXz.T>得到如圖6所示:四、線性回歸擬合對(duì)于多元線性回歸模型y=B+BxH Px+e011pp設(shè)變量Xi,xj…Xp,y盼組觀測(cè)值為(x,x，???兀,y)i=1,2,?…,n?i1 i2 ipi、1

記x=??、1

記x=??1x…xry)rp)121p10x22???…x2p??????,y=y2?，則p=p1?xn2…xnp<y丿njp;x11x21???xn1的估計(jì)值為b=P=(x'x)-1x'y(11.2)在Matlab中，用regress函數(shù)進(jìn)行多元線性回歸分析，應(yīng)用方法如下:語法：b=regress(y,x)[b,bint,r,rint,stats]=regress(y,x)[b,bint,r,rint,stats]=regress(y,x,alpha)b=regress(y,x)得到的p+1維列向量b即為(11.2)式給出的回歸系數(shù)P的估計(jì)值?[b,bint,r,rint,stats]=regress(y,x)給出回歸系數(shù)P的估計(jì)值b，p的95%置信區(qū)間((p+1)x2向量)bint,殘差r以及每個(gè)殘差的95%置信區(qū)間(nx2向量)rint；向量stats給出回歸的R2統(tǒng)計(jì)量和F以及臨界概率p的值.如果卩的置信區(qū)間(bint的第i+1行)不包含0,則在顯著水平為a時(shí)拒絕ip=0的假設(shè)，認(rèn)為變量x是顯著的.ii[b,bint,r,rint,stats]=regress(y,x,alpha)給出了bint和rint的100(l-alpha)%的置信區(qū)間.三次樣條插值函數(shù)的MATLAB程序matlab的splinex=0:10;y=sin(x); %插值點(diǎn)xx=0:.25:10; %繪圖點(diǎn)yy=spline(x,y,xx);plot(x,y,'o',xx,yy)結(jié)果如圖7所示圖7五、非線性擬合非線性擬合可以用以下命令（同樣適用于線形回歸分析）：beta=nlinfit（X,y,fun,betaO）X給定的自變量數(shù)據(jù),Y給定的因變量數(shù)據(jù),fun要擬合的函數(shù)模型（句柄函數(shù)或者內(nèi)聯(lián)函數(shù)形式），beta0函數(shù)模型中系數(shù)估計(jì)初值,beta返回?cái)M合后的系數(shù)x=lsqcurvefit（fun,x0,xdata,ydata）fun要擬合的目標(biāo)函數(shù),x0目標(biāo)函數(shù)中的系數(shù)估計(jì)初值,xdata自變量數(shù)據(jù),ydata函數(shù)值數(shù)據(jù)X擬合返回的系數(shù)（擬合結(jié)果）nlinfit格式：[beta,r，J]=nlinfit（x,y,'model',beta0）Beta估計(jì)出的回歸系數(shù)r殘差JJacobian矩陣x，y輸入數(shù)據(jù)x、y分別為n*m矩陣和n維列向量，對(duì)一元非線性回歸，x為n維列向量。model是事先用m-文件定義的非線性函數(shù)beta0回歸系數(shù)的初值例1已知數(shù)據(jù)：xl二[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ù)a,b,c,d,e。對(duì)回歸模型建立M文件model.m如下：functionyy二myfun(beta,x)x1=x(:,1);%數(shù)據(jù)x2=x(:,2);%數(shù)據(jù)x3=x(:,3);%數(shù)據(jù)yy=beta(1)+beta(2)*x2+beta(3)*x3+beta(4)*(x2."2)+beta(5)*(x3."2);%運(yùn)算出圖主程序如下：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]; %數(shù)據(jù)y=[0.785,0.703,0.583,0.571,0.126];%數(shù)據(jù)beta0=[1,1,1,1,1];%數(shù)據(jù)[beta,r,j]=nlinfit(x,y,@myfun,betaO)%運(yùn)算出圖輸出結(jié)果如圖8所示圖8例題2：混凝土的抗壓強(qiáng)度隨養(yǎng)護(hù)時(shí)間的延長(zhǎng)而增加，現(xiàn)將一批混凝土作成12個(gè)試塊，記錄了養(yǎng)護(hù)日期(日)及抗壓強(qiáng)度y(kg/cm2)的數(shù)據(jù)：養(yǎng)護(hù)時(shí)間：x=[234579121417212856]抗壓強(qiáng)度：y=[35+r42+r47+r53+r59+r65+r68+r73+r76+r82+r86+r99+r]建立非線性回歸模型，對(duì)得到的模型和系數(shù)進(jìn)行檢驗(yàn)。注明：此題中的+r代表加上一個(gè)［-0.5,0.5］之間的隨機(jī)數(shù)模型為：y=a+k1*exp(m*x)+k2*exp(-m*x);Matlab程序：x=[234579121417212856];%數(shù)據(jù)r二rand(l,12)-0.5;y1=[354247535965687376828699];%數(shù)據(jù)y=y1+r;%運(yùn)算路徑myfunc二inline('beta(1)+beta(2)*exp(beta(4)*x)+beta(3)*exp(-beta(4)*x)','beta','x');%運(yùn)算路徑beta=nlinfit(x,y,myfunc,[0.50.50.50.5]);%運(yùn)算方法a二beta(1),k1=beta(2),k2二beta(3),m=beta(4)%數(shù)據(jù)%testthemodel%運(yùn)算方法xx=min(x):max(x);%運(yùn)算方法yy=a+k1*exp(m*xx)+k2*exp(—m*xx);%運(yùn)算方法plot(x,y,'o',xx,yy,'r')結(jié)果如圖9所示:圖9lsqnonlin非線性最小二乘(非線性數(shù)據(jù)擬合)的標(biāo)準(zhǔn)形式為minf(x)=f］(x)2+f2(x)2++fm(x)2+Lx其中：L為常數(shù)在MATLAB5.X中，用函數(shù)leastsq解決這類問題，在6.0版中使用函數(shù)lsqnonlin。

設(shè)F(x)設(shè)F(x)=f1(x)f2(x)fm(x)m則目標(biāo)函數(shù)可表達(dá)為min1||F(x)||2=1工f(x)2x2 2 2ii其中：X為向量，F(xiàn)(x)為函數(shù)向量。函數(shù)lsqnonlin格式x=lsqnonlin(fun,xO)%xO為初始解向量；fun為fi(x),i=l,2,…，m,fun返回向量值F,而不是平方和值，平方和隱含在算法中，fun的定義與前面相同。x=lsqnonlin(fun,xO,lb,ub)%lb、ub定義x的下界和上界：lb<x<ub。x=lsqnonlin(fun,xO,lb,ub,options)%options為扌旨定優(yōu)化參數(shù)，若x沒有界，則lb=[]，ub=[]。[x,resnorm]=lsqnonlin(…)％resnorm二sum(fun(x)."2)，即解x處目標(biāo)函數(shù)值。[x,resnorm,residual]=lsqnonlin(…)％residual二fun(x)，即解x處fun的值。[x,resnorm,residual,exitflag]=lsqnonlin(…)％exitflag為終止迭代條件。[x,resnorm,residual,exitflag,output]二lsqnonlin(…)％output輸出優(yōu)化信息。[x,resnorm,residual,exitflag,output,lambda]=lsqnonlin(…)%lambda為L(zhǎng)agrage乘子。[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(…)%fun在解x處的Jacobian矩陣。例求下面非線性最小二乘問題f(2+2k-ekx1-ekx2)2初始解向量為k=1x0=[0.3,0.4]。解：先建立函數(shù)文件，并保存為myfun.m，由于lsqnonlin中的fun為向量形式而不是平方和形式，因此，myfun函數(shù)應(yīng)由fi(x)建立：fk(x)=2+2k一ekx1一ekx2 k=l,2,…，10functionF=myfun(x)k=1:10;F=2+2*k-exp(k*x(1))-exp(k*x(2));然后調(diào)用優(yōu)化程序：x0=[0.30.4]；[x,resnorm]=lsqnonlin(@myfun,x0)%主運(yùn)算方法結(jié)果為：Optimizationterminatedsuccessfully:NormofthecurrentstepislessthanOPTIONS.TolXx=0.25780.2578resnorm=%求目標(biāo)函數(shù)值lsqcurvefit非線性曲線擬合是已知輸入向量xdata和輸出向量ydata,并且知道輸入與輸出的函數(shù)關(guān)系為ydata二F(x,xdata)，但不知道系數(shù)向量x。今進(jìn)行曲線擬合,求x使得下式成立：min

x(F(x,xdatai)—ydatai)2min

x(F(x,xdatai)—ydatai)2i在MATLAB5.X中，使用函數(shù)curvefit解決這類問題。函數(shù)lsqcurvefit格式x=lsqcurvefit(fun,x0,xdata,ydata)%運(yùn)算方法x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)[x,resnorm]=lsqcurvefit(…)％運(yùn)算方法[x,resnorm,residual]=lsqcurvefit(…)％運(yùn)算方法[x,resnorm,residual,exitflag]=lsqcurvefit(…)％運(yùn)算方法[x,resnorm,residual,exitflag,output]二lsqcurvefit(…)[x,resnorm,residual,exitflag,output,lambda]=lsqcurvefit(…)%運(yùn)算方法[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(…)％運(yùn)算方法參數(shù)說明：xO為初始解向量；xdata,ydata為滿足關(guān)系ydata二F(x,xdata)的數(shù)據(jù);lb、ub為解向量的下界和上界lb<x<ub，若沒有指定界，則lb=[]，ub=[]；options為指定的優(yōu)化參數(shù)；fun為擬合函數(shù)，其定義方式為：x=lsqcurvefit(@myfun,xO,xdata,ydata),其中myfun已定義為 functionF二myfun(x,xdata)F=…％計(jì)算x處擬合函數(shù)值fun的用法與前面相同；resnorm二sum((fun(x,xdata)-ydata)."2),即在x處殘差的平方和；residual二fun(x,xdata)-ydata,即在x處的殘差；exitflag為終止迭代的條件；output為輸出的優(yōu)化信息；lambda為解x處的Lagrange乘子；jacobian為解x處擬合函數(shù)fun的jacobian矩陣。例求解如下最小二乘非線性擬合問題已知輸入向量xdata和輸出向量ydata,且長(zhǎng)度都是n,擬合函數(shù)為ydata(i)=x(1)-xdata(i)2+x(2)-sin(xdata(i))+x(3)-xdata(i)3即目標(biāo)函數(shù)為 xdal)—ydata.)2x 2i=1其中：F(x,xdata)=x(1)-xdata2+x(2)-sin(xdata)+x(3)-xdata3初始解向量為x0=[0.3,0.4,0.1]。解：先建立擬合函數(shù)文件，并保存為myfun.mfunctionF=myfun(x,xdata)

F=x(l)*xdata.2+x(2)*sin(xdata)+x(3)*xdata.3;然后給出數(shù)據(jù)xdata和ydata>>xdata=[3.67.79.34.18.62.81.37.910.05.4];>>ydata=[16.5150.6263.124.7208.59.92.7163.9325.054.3];>>x0=[10,10,10]; %初始估計(jì)值>>[x,resnorm]=lsqcurvefit(@myfun,x0,xdata,ydata)結(jié)果為：Optimizationterminatedsuccessfully:RelativefunctionvaluechangingbylessthanOPTIONS.TolFunx=0.2269 0.3385 0.3021resnorm=6.2950對(duì)于上述這些可化為線性模型的回歸問題，一般先將其化為線性模型，然后再用最小二乘法求出參數(shù)的估計(jì)值，最后再經(jīng)過適當(dāng)?shù)淖儞Q，得到所求回歸曲線。在熟練掌握最小二乘法的情況下，解決上述問題的關(guān)鍵是確定曲線類型和怎樣將其轉(zhuǎn)化為線性模型。如圖10所示圖10亡3圖10亡3整握聲SBW確定曲線類型一般從兩個(gè)方面考慮：一是根據(jù)專業(yè)知識(shí)，從理論上推導(dǎo)或憑經(jīng)驗(yàn)推測(cè)、二是在專業(yè)知識(shí)無能為力的情況下，通過繪制和觀測(cè)散點(diǎn)圖確定曲線大體類型。六、Subplot問題程序：x1=[1 1 1 1 1 2 2 2 2 3 333 4 4 4 4 5 5 5 6 6 6 6 7 8 8 8 8 1010 101011 1112 1213 1314 1516 1616 1720]';%數(shù)據(jù)x2=[1 0 1 0 0 1 0 0 0 0 111 0 1 0 0 0 0 1 0 1 0 1 1 0 1 1 0 001 11001 0 1 0 1 1 0 0 0 0]';%數(shù)據(jù)x5二[1,3,3,2,3,2,2,1,3,2,1,2,3,1,3,3,2,2,3,1,1,3,2,2,1,2,1,3,1,1,2,3,2,2,1,2,3,1,2,2,3,2,2,1,2,1]';%數(shù)據(jù)x3=[1000000100100100000110001010110000100100000101]';%數(shù)據(jù)x4=[0001 011 001010000 110 000 110 100001011010011 011 010]';%數(shù)據(jù)y=[13876116081870111283117672087211772 10535 121951231314975213711980011417 20263 132311288413245 13677 15965 12366 2135213839228841697814803 17404 221841354814467 15942 23174 23780 2541014861168822417015990 26330 179492568527837 18838 17483 19207 19346];%數(shù)據(jù)size(x5)x=[ones(46,1),x1,x2,x3,x4][b,bint,r,rint,stats]二regress(y,x,0.05)%運(yùn)算方法plot(x1,r,'r+')holdonfori=1:1:46if(x2(i,1)==0&x5(i,1)==1)k(i)=1;%交替運(yùn)算elseif(x2(i,1)==1&x5(i,1)==1)k(i)=2;%交替運(yùn)算elseif(x2(i,1)==0&x5(i,1)==2)k(i)=3;%交替運(yùn)算elseif(x2(i,1)==1&x5(i,1)==2)k(i)=4;%交替運(yùn)算elseif(x2(i,1)==0&x5(i,1)==3)k(i)=5;%交替運(yùn)算elsek(i)=6;endendk;plot(k,r,'g+')輸出如圖11所示圖11第二個(gè)TOC\o"1-5"\h\zx1=[11 1 1 1 2 2 2 2 3 333 4 4 4 4 5 5 5 6 6 6 6 7 8 8 8 8 1010 1010111112 1213 1314 1516 16 16 1720]';%數(shù)據(jù)x2=[10 1 0 0 1 0 0 0 0 111 0 1 0 0 0 0 1 0 1 0 1 1 0 1 1 0 001 110010 1 0 1 1 0 0 0 0]';%數(shù)據(jù)x5二[1,3,3,2,3,2,2,1,3,2,1,2,3,1,3,3,2,2,3,1,1,3,2,2,1,2,1,3,1,1,2,3,2,2,1,2,3,1,2,2,3,2,2,1,2,1]';%數(shù)據(jù)x3=[10 0 0 0 0 0 1 0 01001 0 0 0 0 0 1 1 0 0 0 1 0 1 0 1 1 0 0 0 01001 0 0 0 0 0 1 0 1]';%數(shù)據(jù)x4=[00 0 1 0 1 1 0 0 10100 0 0 1 1 0 0 0 0 1 1 0 1 0 0 0 0 1 0 1 10100 1 1 0 1 1 0 1 0]';%數(shù)據(jù)y二[1387611608187011128311