matlab在統(tǒng)計(jì)中的應(yīng)用

1、數(shù)理統(tǒng)計(jì)建模Matlab在統(tǒng)計(jì)中的應(yīng)用,山西財(cái)經(jīng)大學(xué)應(yīng)用數(shù)學(xué)學(xué)院 高崇山,一、 概率分布及有關(guān)函數(shù),調(diào)用格式為：分布命令符功能命令符(x,參數(shù))[m,v]=分布stat(x,參數(shù))rand產(chǎn)生[0,1]上的隨機(jī)數(shù)，randn產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)。,y=normpdf(x,mu,sigma) 或y=pdf(‘norm’,x,mu,sigma) % 正態(tài)分布N(mu,sigma2)在x處的概率密度；y=

2、normcdf(x,mu,sigma) 或y=cdf(‘norm’,x,mu,sigma) % 正態(tài)分布N(mu,sigma2)在x處的分布函數(shù)；y=norminv(alpha,mu,sigma) %正態(tài)分布N(mu,sigma2)在對(duì)應(yīng)于alfa的分位數(shù)。即 [m,v]=normstat(mu,sigma) %正態(tài)分布N(mu,sigma2)的期望和方差；y=exprnd(lamda,[m,n])或random(‘

3、exp’,lamda,[m,n]) %產(chǎn)生一個(gè)m × n的服從參數(shù)為lamda的指數(shù)分布的隨機(jī)矩陣,二、 描述性統(tǒng)計(jì) 描述性統(tǒng)計(jì)就是搜集、整理、加工和分析統(tǒng)計(jì)數(shù)據(jù)，使之系統(tǒng)化，以顯示出數(shù)據(jù)的趨勢(shì)、特征和數(shù)量關(guān)系。,2.1 樣本均值mean和中值median 它們都是樣本數(shù)據(jù)在數(shù)據(jù)分布線(xiàn)上中心位置的度量.,A=[1 2 4 4;3 4 6 6;5 6 8 8;5 6 8 8];mean

4、(A) % 計(jì)算矩陣每列的均值，相當(dāng)于mean(A,1)mean(A,2) % 計(jì)算矩陣每行的均值median(A) % 計(jì)算矩陣每列的中值（中位數(shù)），相當(dāng)于median(A,1)median(A,2) % 計(jì)算矩陣每行的中值（中位數(shù)）,2.2方差var、標(biāo)準(zhǔn)差std、極差range和協(xié)方差cov 它們都是描述樣本中的數(shù)據(jù)偏離其中心值的程度,X=rand(4,5);

5、std(X) % 計(jì)算矩陣X每列的標(biāo)準(zhǔn)差var(X) %計(jì)算矩陣X每列的方差range(X) %計(jì)算矩陣X每列的極差cov(X) %計(jì)算協(xié)方差var(X)=diag(cov(X))’std(X)=sqrt(diag(cov(X)))’,X若為向量，cov(X)=var(X)；若X為矩陣，X的每一列表示一個(gè)變量而行元素為觀察值。,對(duì)于二維隨機(jī)向量(X,Y)，x為X的觀察值，y為Y

6、的觀察值(x,y為同維向量)，則有：cov(x,y)=cov([x,y]),2.3 百分位數(shù)及其圖形描述 百分位數(shù)(percentile)是把數(shù)據(jù)按從小到大的順序排列后，位于p%位置的值稱(chēng)為第p百分位數(shù)。第25百分位數(shù)由叫做四分之一分位數(shù)(下四分位數(shù))，75百分位數(shù)由叫做四分之三分位數(shù)(上四分位數(shù))， 第50百分位數(shù)就是median中數(shù)。最小值是第0百分位數(shù)，最大值是第100百分位數(shù)。 百分

7、位數(shù)是用于反映樣本數(shù)據(jù)形態(tài)信息的數(shù)據(jù)統(tǒng)計(jì)量，它也可以刻劃數(shù)據(jù)的位置和散布特征。 Y = prctile(X,p) 返回樣本X中大于p%(0<p<100)的值。如果X是向量，則返回X中p百分位數(shù)，若X為矩陣，則返回一個(gè)關(guān)于每列元素的p百分位數(shù)行向量。注意：p也可以是一個(gè)向量，此時(shí)，返回一組百分位數(shù)。,eg. x =100*rand(1,10) ,y=prctile(x,0:10:100),s

8、ubplot(1,2,1),boxplot(x),subplot(1,2,2),bar(x),x=[61.5432 79.1937 92.1813 73.8207 17.6266 40.5706 93.5470 91.6904 41.0270 89.3650],y=[17.6266 29.0986 40.7988 51.2851 67.6820 76.5072 84.279

9、3 90.5277 91.9359 92.8641 93.5470],2.4 相關(guān)系數(shù) 相關(guān)系數(shù)反映兩個(gè)隨即變量之間線(xiàn)性相依程度的變量。R=corrcoef(X)R=corrcoef(x,y)含義同協(xié)方差。同協(xié)方差之間有如下關(guān)系,2.5 樣本峰度和偏度偏度描述的是分布的對(duì)稱(chēng)性，其定義為： 當(dāng)f >0時(shí)，表

10、示數(shù)據(jù)在均值右邊的比左邊的多；f <0正好相反；f 接近于0，則表示分布是對(duì)稱(chēng)的。峰度描述的是分布曲線(xiàn)的陡緩程度，定義為： 它是以正態(tài)分布為標(biāo)準(zhǔn)，比較兩側(cè)極端數(shù)據(jù)分布的情況的指標(biāo)。g較大，則表示樣本中有許多遠(yuǎn)離均值的數(shù)據(jù)。上述公式中，s是樣本標(biāo)準(zhǔn)差。f=skewness(X)g=kurtosis(X),三、 參數(shù)估計(jì),1）fit函數(shù)的調(diào)用方法類(lèi)似，以正態(tài)分布說(shuō)明之。格式為：[mu,sigma,muc

11、i,sigmaci]=normfit(x,alpha)說(shuō)明：x是樣本(矩陣或向量)；alpha是顯著性水平(默認(rèn)值為0.05)；mu是總體均值的點(diǎn)估計(jì)值；sigma是總體方差的點(diǎn)估計(jì)值；muci是總體均值的區(qū)間估計(jì)；sigmaci是總體方差的區(qū)間估計(jì)。2）mle的調(diào)用格式[phat,pci]=mle(‘dist’,data,alpha,pl)說(shuō)明：dist是所給的分布名(如:norm,exp…)；data是樣本數(shù)據(jù)；alpha為

12、可選項(xiàng)，表示顯著性水平；pl僅用于二項(xiàng)分布，表示試驗(yàn)的次數(shù)；phat為返回的點(diǎn)估計(jì)值；pci為返回相應(yīng)置信區(qū)間。,eg. x=normrnd(2,4,100,1); [mu,sigma,muci,sigmaci]=normfit(x) [mu,sigma,muci,sigmaci]=normfit(x,0.1) [phat,pci]=mle('norm',x) [phat,

13、pci]=mle('norm',x,0.1),四、 假設(shè)經(jīng)驗(yàn)4.1 單個(gè)樣本的t檢驗(yàn)功能：進(jìn)行樣本均值的t檢驗(yàn)格式：h=ttest(x,m)；% 在0.05的顯著性水平下進(jìn)行t檢驗(yàn)，以確定在標(biāo)準(zhǔn)差未知的情況下取自正態(tài)分布的均值是否為m，若輸出h=0，則接受零假設(shè)，h=1則否定零假設(shè)；h=ttest(x,m,alpha)； % alpha 為給定的顯著性水平。[h,sig,ci]=ttest(x,m,alp

14、ha,tail) % 若原假設(shè)為μ=μ0，則取tail=0(可省略)；若原假設(shè)為μ>μ0，則取tail=1；若原假設(shè)為μalpha時(shí)不能否定零假設(shè)，一般sig越大零假設(shè)越可信。ci為均值真值的1-alpha置信區(qū)間。,4.2 單個(gè)樣本的z檢驗(yàn)功能：在給定方差的條件下進(jìn)行z檢驗(yàn)格式：h=ztest(x,m,sigma) % sigma 正態(tài)總體的標(biāo)準(zhǔn)差，alpha = 0.05h=ztest(x,m,sigma,a

15、lpha)[h,sig,ci]=ztest(x,m,sigma,alpha,tail)4.3 兩個(gè)樣本的t檢驗(yàn)功能：兩個(gè)服從正態(tài)總體樣本均值差異的t檢驗(yàn)(σ12 = σ22均未知)格式： [h,significance,ci] = ttest2(x,y) % 默認(rèn)alpha =0.05 [h,significance,ci] = ttest2(x,y,alpha) [h,significance,ci] = t

16、test2(x,y,alpha,tail),五、 統(tǒng)計(jì)繪圖5.1 box圖boxplot(X,notch,’sym’,vert,whis)X：樣本數(shù)據(jù)；notch=1有切口, notch=0無(wú)切口,(默認(rèn)notch=0)；‘sym’野值標(biāo)記符號(hào),默認(rèn)‘+’；vert=0,box圖是水平放置,vert=1是垂直放置(默認(rèn))。whis 定義虛線(xiàn)的長(zhǎng)度，一般用缺省值.,eg:x1=normrnd(4,1,200,1);

17、x2=normrnd(8,1,200,1);x3=normrnd(6,2,200,1);x=[x1,x2,x3];boxplot(x,1)圖見(jiàn)下頁(yè),說(shuō)明：1.盒子的上下兩條線(xiàn)分別為樣本的75%和25%分位線(xiàn)，中間為樣本中位數(shù)；2.虛線(xiàn)表示樣本的其余部分，位于盒子的上下兩側(cè)；3.‘+’表示野值（奇異值），位于虛線(xiàn)的上方和下方；4.‘切口’表示樣本中位數(shù)的置信區(qū)間。默認(rèn)狀態(tài)下無(wú)切口。,5.2 正態(tài)概率圖正態(tài)概率圖用于判斷

18、樣本數(shù)據(jù)是否服從正態(tài)分布。格式：normplot(X)X：數(shù)據(jù).若X為矩陣，則為X的每列顯示一條線(xiàn)。圖形以符號(hào)‘+’顯示樣本數(shù)據(jù)。如果數(shù)據(jù)服從正態(tài)分布，則圖形呈現(xiàn)直線(xiàn)，否則會(huì)表現(xiàn)不同程度的曲線(xiàn)。,x1=normrnd(4,1,200,1);x2=normrnd(8,1,200,1);x3=normrnd(6,2,200,1);x=[x1,x2,x3];normplot(x,1),5.3 分位數(shù)—分位數(shù)圖分位數(shù)—分位數(shù)圖用于比較

19、兩個(gè)樣本的分布.格式:qqplot(X,Y,pvec)其中，X,Y分別是兩個(gè)樣本的數(shù)據(jù)。如果兩個(gè)樣本來(lái)自同一分布，則繪制的曲線(xiàn)為直線(xiàn)。若X,Y為矩陣，則為他們每一列顯示一條直線(xiàn)。圖形以符號(hào)‘+’顯示樣本數(shù)據(jù)。參數(shù)pvec是可選項(xiàng)，用于規(guī)定分位數(shù)。,x1=normrnd(4,1,200,1);x2=normrnd(0,1,200,1); qqplot(x1,x2),六、 分布檢驗(yàn)6.1 Jarque-bera檢驗(yàn)該檢驗(yàn)評(píng)價(jià)

20、X服從未知均值和方差的正態(tài)分布的假設(shè)是否成立。該檢驗(yàn)基于X的樣本偏度和峰度。對(duì)于正態(tài)分布數(shù)據(jù)，偏度接近于0，峰度接近于3。 Jarque-bera檢驗(yàn)就是確定樣本偏度、峰度是否與它們的期望值相差較遠(yuǎn)。功能：測(cè)試數(shù)據(jù)對(duì)正態(tài)分布的擬合程度。格式:h = jbtest(X) % 當(dāng)h=1,拒絕X服從正態(tài)分布;否則h=0。(默認(rèn)alpha=0.05)h = jbtest(X,alpha) [h,P,jbs

21、tat,cv] = jbtest(X,alpha) % P為檢驗(yàn)的p值，jbstat為檢驗(yàn)的統(tǒng)計(jì)量，cv為確定是否拒絕零假設(shè)的臨界值。當(dāng)jbstat<cv時(shí)，同樣接受零假設(shè)。注意：該檢驗(yàn)不能用于小樣本的檢驗(yàn)，只能用于大樣本。對(duì)于小樣本，用lillietest檢驗(yàn)較合適。,6.2 Lilliefors檢驗(yàn)該檢驗(yàn)評(píng)價(jià)X服從未知均值和方差的正態(tài)分布的假設(shè)是否成立，對(duì)應(yīng)的備擇假設(shè)為X不服從正態(tài)分布。本檢驗(yàn)比較X的經(jīng)

22、驗(yàn)分布與具有相同均值和方差的正態(tài)分布。格式：H=lillietest(X) % 若H=0，則接受X服從正態(tài)分布；否則，H=1。(默認(rèn)alpha=0.05)H=lillietest(X,alpha)[H,P,LSTAT,CV]= lillietest(X,alpha) % P為檢驗(yàn)的p值，通過(guò)在一系列由Lilliefors創(chuàng)建的表中進(jìn)行插值得到；LSTAT為檢驗(yàn)統(tǒng)計(jì)量的值；CV為確定是否拒絕零假設(shè)的臨界值。如果LST

23、AT的值位于Lilliefors表之外，則P返回NaN，但H顯示是否拒絕假設(shè)。當(dāng)LSTAT>CV時(shí)，同樣拒絕零假設(shè)。,七、 回歸分析7.1多元線(xiàn)性回歸分析數(shù)學(xué)模型一元回歸模型為：y = β0+ β1x + ε 其中ε服從N(0,1)多元回歸模型為： y = β0+ β1x1+ β2x2+…+ βmxm+ ε 其中ε服從N(0,σ2)

24、 回歸問(wèn)題就是求出xi的系數(shù)βi，并求出誤差σ2的估計(jì)，回歸系數(shù)β的區(qū)間估計(jì)和假設(shè)檢驗(yàn)，模型的有效性檢驗(yàn)及對(duì)給定的x做出y的預(yù)測(cè)。 預(yù)測(cè)分點(diǎn)預(yù)測(cè)和區(qū)間預(yù)測(cè)：其中點(diǎn)預(yù)測(cè)將x代入模型中即可，區(qū)間預(yù)測(cè)需要編一個(gè)小程序。,命令為：b=regress(y,x)[b,bint,r,rint,s]=regress(y,x,alpha)說(shuō)明： 輸入y(因變量，列向量)；x(第一列全為1，第二

25、列為x1的觀察值，第三列為x2的觀察值，…)；alpha為顯著性水平α(默認(rèn)值為0.05)； 輸出b為(β0, β2,… β2, βm)的估計(jì)值；bint是(β0, β2,… β2, βm)的置信區(qū)間；r是殘差(觀察值與預(yù)測(cè)值之差，為列向量，主要用于探測(cè)模型假設(shè)的合理性)，rint是殘差的置信區(qū)間；s包含3個(gè)統(tǒng)計(jì)量：第一個(gè)是決定系數(shù)R2(其值越大，說(shuō)明自變量對(duì)因變量的所起的作用也越大，但無(wú)明確界限說(shuō)明模型是否有效)，第二個(gè)是F

26、值，第三個(gè)是F(1,N-2)分布大于F值的概率p，p<alpha時(shí)，回歸模型有效。,y=[144215138145162142170124158154 162150140110128130135114116124 136142120120160158144130125175];x=[39474547654667426756 645659344248

27、45182019 36503921445363292569];n=length(y);X=[ones(n,1) x'];[b,bint,r,rint,s]=regress(y',X);b,bint,s,s2=sum(r.^2)/(n-2)rcoplot(r,rint)pausey=[y(1) y(3:30)];x=[x(1) x(3:30)];n=length(y);X=[

28、ones(n,1) x'];[b,bint,r,rint,s]=regress(y',X);,b,bint,s,s2=sum(r.^2)/(n-2)rcoplot(r,rint)pausey0=b(1)+b(2)*50; % 預(yù)測(cè)y(x=50)xb=mean(x);sxx=sum((x-xb).^2);a=sqrt((50-xb)^2/sxx+1/n+1);t=tinv(0.9

29、75,n-2);d=t*a*sqrt(s2);y1=y0-d;y2=y0+d; % 預(yù)測(cè)y(x=50)區(qū)間（t分布）[y0 y1 y2]d1=norminv(0.975)*sqrt(s2);y3=y0-d1;y4=y0+d1;[y0 y3 y4] % 預(yù)測(cè)y(x=50)區(qū)間（N分布）,7.2 多項(xiàng)式回歸1)多項(xiàng)式曲線(xiàn)擬合 多項(xiàng)式回歸的模型為：p(x)=p1

30、xn+p2xn-2+…+pnx+pn+1 格式為：[p,s]=polyfit(x,y,n)其中：n為擬合次數(shù)；x,y分別是自變量和因變量；s是一個(gè)矩陣，用于polyval函數(shù)，可進(jìn)行預(yù)測(cè)的誤差估計(jì)；p為系數(shù)向量( p1,p2,…,pn,pn+1 )的估計(jì)值。2）多項(xiàng)式預(yù)測(cè)和置信區(qū)間的評(píng)估 [y,delta]=polyconf(p,x,s,alpha)其中：p,s是擬合輸出的結(jié)果，x是要預(yù)測(cè)的點(diǎn)，alpha是置信度，輸出

31、的是1-alpha的置信區(qū)間y±delta.說(shuō)明：命令polytool(x,y,n,alpha)作用類(lèi)似于polyfit.他是一個(gè)交互式畫(huà)面。,eg. y=[1035 624 1084 1052 1015 1066 704 960 990 1050 839 1030 985 855]

32、;x=[6.0000 2.5000 7.5000 8.5000 10.0000 7.0000 3.0000 11.5000 5.5000 6.5000 4.0000 9.0000 11.0000 12.5000];plot(x,y,'+'),pausex2=x.^2;X=[ones(14,1) x' x2'];[b,bi,r,

33、ri,s]=regress(y',X);b,bi,s,pausexx=2:.1:13;yy=b(1)+b(2)*xx+b(3)*xx.^2;plot(x,y,'+',xx,yy),grid,pausea=polyfit(x,y,2) % 注意與 regress 的區(qū)別polytool(x,y,2),7.3 多元二項(xiàng)式回歸模型： y = b0 + b1x1+b2x2+b3

34、x3+… (線(xiàn)性項(xiàng)) + b11x12+b22x22+b33x32+… (二次項(xiàng)) +b12x1x2+b13x1x3+…. (交叉項(xiàng))格式：rstool(x,y) %用于顯示模型的交互式預(yù)測(cè)圖，包括預(yù)測(cè)的95%置信區(qū)間。rstool(x,y,’model’,alpha)

35、 alpha為可選項(xiàng)。% model=‘interaction’表示回歸模型包括常數(shù)項(xiàng)、線(xiàn)性項(xiàng)、交叉項(xiàng)；% model=‘quadratic’表示回歸模型包括常數(shù)項(xiàng)、線(xiàn)性項(xiàng)、交叉項(xiàng)、二次項(xiàng)；% model=‘purequadratic’表示回歸模型包括常數(shù)項(xiàng)、線(xiàn)性項(xiàng)、二次項(xiàng)；rstool(x,y,’model’,alpha,’xname’,’yname’)% 可以在x軸和y軸上分別標(biāo)上相應(yīng)的名稱(chēng)。,eg.n=10;x

36、1=[120 140 190 130 155 175 125 145 180 150];x2=[100 110 90 150 210 150 250 270 300 250];y =[102 100 120 77 46 93 26 69 65 85];x=[ones(n,1) x1' x2'];[b,bi,r,ri,s]=regress(y',x);s2=sum(r.^2)/(n-2);