Matlab在統(tǒng)計中的應(yīng)用

數(shù)理統(tǒng)計建模

Matlab在統(tǒng)計中的應(yīng)用山西財經(jīng)大學(xué)應(yīng)用數(shù)學(xué)學(xué)院高崇山一、概率分布及有關(guān)函數(shù)字符分布unif均勻分布exp指數(shù)分布norm正態(tài)分布chi2chi方分布tt分布ff分布bino二項分布poisspoisson分布字符功能調(diào)用格式pdf概率密度namepdf(x,參數(shù)）cdf分布函數(shù)inv逆概率分布stat均值與方差ran隨機數(shù)生成調(diào)用格式為：分布命令符功能命令符(x,參數(shù))[m,v]=分布stat(x,參數(shù))rand產(chǎn)生[0,1]上的隨機數(shù)，randn產(chǎn)生標(biāo)準(zhǔn)正態(tài)分布隨機數(shù)。y=normpdf(x,mu,sigma)或y=pdf(‘norm’,x,mu,sigma)%正態(tài)分布N(mu,sigma2)在x處的概率密度；y=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)在對應(yīng)于alfa的分位數(shù)。即[m,v]=normstat(mu,sigma)%正態(tài)分布N(mu,sigma2)的期望和方差；y=exprnd(lamda,[m,n])或random(‘exp’,lamda,[m,n])%產(chǎn)生一個m×n的服從參數(shù)為lamda的指數(shù)分布的隨機矩陣二、描述性統(tǒng)計描述性統(tǒng)計就是搜集、整理、加工和分析統(tǒng)計數(shù)據(jù)，使之系統(tǒng)化，以顯示出數(shù)據(jù)的趨勢、特征和數(shù)量關(guān)系。函數(shù)名描述函數(shù)名描述max求向量或矩陣列的最大值sort升序排列min求向量或矩陣列的最小值sum求向量或矩陣列的和mean求向量或矩陣列的平均值cumsum累計求和median求向量或矩陣列的中間值cov求協(xié)方差std求標(biāo)準(zhǔn)差corrcoef求相關(guān)系數(shù)var求方差kurtosis計算樣本峰度rang樣本極差skewness計算樣本偏度2.1樣本均值mean和中值median它們都是樣本數(shù)據(jù)在數(shù)據(jù)分布線上中心位置的度量.A=[1244;3466;5688;5688];mean(A)%計算矩陣每列的均值，相當(dāng)于mean(A,1)mean(A,2)%計算矩陣每行的均值median(A)%計算矩陣每列的中值〔中位數(shù)〕，相當(dāng)于median(A,1)median(A,2)%計算矩陣每行的中值〔中位數(shù)〕2.2方差var、標(biāo)準(zhǔn)差std、極差range和協(xié)方差cov它們都是描述樣本中的數(shù)據(jù)偏離其中心值的程度X=rand(4,5);std(X)%計算矩陣X每列的標(biāo)準(zhǔn)差var(X)%計算矩陣X每列的方差range(X)%計算矩陣X每列的極差cov(X)%計算協(xié)方差var(X)=diag(cov(X))’std(X)=sqrt(diag(cov(X)))’X假設(shè)為向量，cov(X)=var(X)；假設(shè)X為矩陣，X的每一列表示一個變量而行元素為觀察值。對于二維隨機向量(X,Y)，x為X的觀察值，y為Y的觀察值(x,y為同維向量)，那么有：cov(x,y)=cov([x,y])2.3百分位數(shù)及其圖形描述百分位數(shù)(percentile)是把數(shù)據(jù)按從小到大的順序排列后，位于p%位置的值稱為第p百分位數(shù)。第25百分位數(shù)由叫做四分之一分位數(shù)(下四分位數(shù))，75百分位數(shù)由叫做四分之三分位數(shù)(上四分位數(shù))，第50百分位數(shù)就是median中數(shù)。最小值是第0百分位數(shù)，最大值是第100百分位數(shù)。百分位數(shù)是用于反映樣本數(shù)據(jù)形態(tài)信息的數(shù)據(jù)統(tǒng)計量，它也可以刻劃數(shù)據(jù)的位置和散布特征。Y=prctile(X,p)返回樣本X中大于p%(0<p<100)的值。如果X是向量，那么返回X中p百分位數(shù)，假設(shè)X為矩陣，那么返回一個關(guān)于每列元素的p百分位數(shù)行向量。注意：p也可以是一個向量，此時，返回一組百分位數(shù)。eg.x=100*rand(1,10),y=prctile(x,0:10:100),subplot(1,2,1),boxplot(x),subplot(1,2,2),bar(x)x=[61.543279.193792.181373.820717.626640.570693.547091.690441.027089.3650]y=[17.626629.098640.798851.285167.682076.507284.279390.527791.935992.864193.5470]2.4相關(guān)系數(shù)

相關(guān)系數(shù)反映兩個隨即變量之間線性相依程度的變量。R=corrcoef(X)R=corrcoef(x,y)含義同協(xié)方差。同協(xié)方差之間有如下關(guān)系2.5樣本峰度和偏度偏度描述的是分布的對稱性，其定義為：當(dāng)f>0時，表示數(shù)據(jù)在均值右邊的比左邊的多；f<0正好相反；f接近于0，那么表示分布是對稱的。峰度描述的是分布曲線的陡緩程度，定義為：它是以正態(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ù)估計函數(shù)名功能binofit二項分布參數(shù)估計expfit指數(shù)分布的參數(shù)估計normfit正態(tài)分布的參數(shù)估計poissfitpoisson分布的參數(shù)估計unifit均勻分布的參數(shù)估計mle最大似然估計……1〕fit函數(shù)的調(diào)用方法類似，以正態(tài)分布說明之。格式為：[mu,sigma,muci,sigmaci]=normfit(x,alpha)說明：x是樣本(矩陣或向量)；alpha是顯著性水平(默認(rèn)值為0.05)；mu是總體均值的點估計值；sigma是總體方差的點估計值；muci是總體均值的區(qū)間估計；sigmaci是總體方差的區(qū)間估計。2〕mle的調(diào)用格式[phat,pci]=mle(‘dist’,data,alpha,pl)說明：dist是所給的分布名(如:norm,exp…)；data是樣本數(shù)據(jù)；alpha為可選項，表示顯著性水平；pl僅用于二項分布，表示試驗的次數(shù)；phat為返回的點估計值；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,pci]=mle('norm',x,0.1)四、假設(shè)經(jīng)驗4.1單個樣本的t檢驗功能：進(jìn)行樣本均值的t檢驗格式：h=ttest(x,m)；%在0.05的顯著性水平下進(jìn)行t檢驗，以確定在標(biāo)準(zhǔn)差未知的情況下取自正態(tài)分布的均值是否為m，假設(shè)輸出h=0，那么接受零假設(shè)，h=1那么否認(rèn)零假設(shè)；h=ttest(x,m,alpha)；%alpha為給定的顯著性水平。[h,sig,ci]=ttest(x,m,alpha,tail)%假設(shè)原假設(shè)為μ=μ0，那么取tail=0(可省略)；假設(shè)原假設(shè)為μ>μ0，那么取tail=1；假設(shè)原假設(shè)為μ<μ0，那么取tail=-1。sig表示在零假設(shè)下，樣本均值出現(xiàn)的概率，當(dāng)sig>alpha時不能否認(rèn)零假設(shè)，一般sig越大零假設(shè)越可信。ci為均值真值的1-alpha置信區(qū)間。4.2單個樣本的z檢驗功能：在給定方差的條件下進(jìn)行z檢驗格式：h=ztest(x,m,sigma)%sigma正態(tài)總體的標(biāo)準(zhǔn)差，alpha=0.05h=ztest(x,m,sigma,alpha)[h,sig,ci]=ztest(x,m,sigma,alpha,tail)4.3兩個樣本的t檢驗功能：兩個服從正態(tài)總體樣本均值差異的t檢驗(σ12=σ22均未知)格式：[h,significance,ci]=ttest2(x,y)%默認(rèn)alpha=0.05[h,significance,ci]=ttest2(x,y,alpha)[h,significance,ci]=ttest2(x,y,alpha,tail)五、統(tǒng)計繪圖5.1box圖boxplot(X,notch,’sym’,vert,whis)X：樣本數(shù)據(jù)；notch=1有切口,notch=0無切口,(默認(rèn)notch=0)；‘sym’野值標(biāo)記符號,默認(rèn)‘+’；vert=0,box圖是水平放置,vert=1是垂直放置(默認(rèn))。whis定義虛線的長度，一般用缺省值.eg:x1=normrnd(4,1,200,1);x2=normrnd(8,1,200,1);x3=normrnd(6,2,200,1);x=[x1,x2,x3];boxplot(x,1)圖見下頁說明：1.盒子的上下兩條線分別為樣本的75%和25%分位線，中間為樣本中位數(shù)；2.虛線表示樣本的其余局部，位于盒子的上下兩側(cè)；3.‘+’表示野值〔奇異值〕，位于虛線的上方和下方；4.‘切口’表示樣本中位數(shù)的置信區(qū)間。默認(rèn)狀態(tài)下無切口。5.2正態(tài)概率圖正態(tài)概率圖用于判斷樣本數(shù)據(jù)是否服從正態(tài)分布。格式：normplot(X)X：數(shù)據(jù).假設(shè)X為矩陣，那么為X的每列顯示一條線。圖形以符號‘+’顯示樣本數(shù)據(jù)。如果數(shù)據(jù)服從正態(tài)分布，那么圖形呈現(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ù)圖用于比較兩個樣本的分布.格式:qqplot(X,Y,pvec)其中，X,Y分別是兩個樣本的數(shù)據(jù)。如果兩個樣本來自同一分布，那么繪制的曲線為直線。假設(shè)X,Y為矩陣，那么為他們每一列顯示一條直線。圖形以符號‘+’顯示樣本數(shù)據(jù)。參數(shù)pvec是可選項，用于規(guī)定分位數(shù)。x1=normrnd(4,1,200,1);x2=normrnd(0,1,200,1);qqplot(x1,x2)六、分布檢驗6.1Jarque-bera檢驗該檢驗評價X服從未知均值和方差的正態(tài)分布的假設(shè)是否成立。該檢驗基于X的樣本偏度和峰度。對于正態(tài)分布數(shù)據(jù)，偏度接近于0，峰度接近于3。Jarque-bera檢驗就是確定樣本偏度、峰度是否與它們的期望值相差較遠(yuǎn)。功能：測試數(shù)據(jù)對正態(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,jbstat,cv]=jbtest(X,alpha)%P為檢驗的p值，jbstat為檢驗的統(tǒng)計量，cv為確定是否拒絕零假設(shè)的臨界值。當(dāng)jbstat<cv時，同樣接受零假設(shè)。注意：該檢驗不能用于小樣本的檢驗，只能用于大樣本。對于小樣本，用lillietest檢驗較適宜。6.2Lilliefors檢驗該檢驗評價X服從未知均值和方差的正態(tài)分布的假設(shè)是否成立，對應(yīng)的備擇假設(shè)為X不服從正態(tài)分布。本檢驗比較X的經(jīng)驗分布與具有相同均值和方差的正態(tài)分布。格式：H=lillietest(X)%假設(shè)H=0，那么接受X服從正態(tài)分布；否那么，H=1。(默認(rèn)alpha=0.05)H=lillietest(X,alpha)[H,P,LSTAT,CV]=lillietest(X,alpha)%P為檢驗的p值，通過在一系列由Lilliefors創(chuàng)立的表中進(jìn)行插值得到；LSTAT為檢驗統(tǒng)計量的值；CV為確定是否拒絕零假設(shè)的臨界值。如果LSTAT的值位于Lilliefors表之外，那么P返回NaN，但H顯示是否拒絕假設(shè)。當(dāng)LSTAT>CV時，同樣拒絕零假設(shè)。七、回歸分析7.1多元線性回歸分析數(shù)學(xué)模型一元回歸模型為：y=β0+β1x+ε其中ε服從N(0,1)多元回歸模型為：

y=β0+β1x1+β2x2+…+βmxm+ε

其中ε服從N(0,σ2)

回歸問題就是求出xi的系數(shù)βi，并求出誤差σ2的估計，回歸系數(shù)β的區(qū)間估計和假設(shè)檢驗，模型的有效性檢驗及對給定的x做出y的預(yù)測。預(yù)測分點預(yù)測和區(qū)間預(yù)測：其中點預(yù)測將x代入模型中即可，區(qū)間預(yù)測需要編一個小程序。命令為：b=regress(y,x)[b,bint,r,rint,s]=regress(y,x,alpha)說明：

輸入y(因變量，列向量)；x(第一列全為1，第二列為x1的觀察值，第三列為x2的觀察值，…)；alpha為顯著性水平α(默認(rèn)值為0.05)；輸出b為(β0,β2,…β2,βm)的估計值；bint是(β0,β2,…β2,βm)的置信區(qū)間；r是殘差(觀察值與預(yù)測值之差，為列向量，主要用于探測模型假設(shè)的合理性)，rint是殘差的置信區(qū)間；s包含3個統(tǒng)計量：第一個是決定系數(shù)R2(其值越大，說明自變量對因變量的所起的作用也越大，但無明確界限說明模型是否有效)，第二個是F值，第三個是F(1,N-2)分布大于F值的概率p，p<alpha時，回歸模型有效。y=[144 215 138 145 162 142 170 124 158 154162 150 140 110 128 130 135 114 116 124136 142 120 120 160 158 144 130 125 175];x=[39 47 45 47 65 46 67 42 67 5664 56 59 34 42 48 45 18 20 1936 50 39 21 44 53 63 29 25 69];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=[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ù)測y(x=50)xb=mean(x);sxx=sum((x-xb).^2);a=sqrt((50-xb)^2/sxx+1/n+1);t=tinv(0.975,n-2);d=t*a*sqrt(s2);y1=y0-d;y2=y0+d;%預(yù)測y(x=50)區(qū)間〔t分布〕[y0y1y2]d1=norminv(0.975)*sqrt(s2);y3=y0-d1;y4=y0+d1;[y0y3y4]%預(yù)測y(x=50)區(qū)間〔N分布〕7.2多項式回歸1)多項式曲線擬合多項式回歸的模型為：p(x)=p1xn+p2xn-2+…+pnx+pn+1格式為：[p,s]=polyfit(x,y,n)其中：n為擬合次數(shù)；x,y分別是自變量和因變量；s是一個矩陣，用于polyval函數(shù)，可進(jìn)行預(yù)測的誤差估計；p為系數(shù)向量(p1,p2,…,pn,pn+1)的估計值。2〕多項式預(yù)測和置信區(qū)間的評估[y,delta]=polyconf(p,x,s,alpha)其中：p,s是擬合輸出的結(jié)果，x是要預(yù)測的點，alpha是置信度，輸出的是1-alpha的置信區(qū)間y±delta.說明：命令polytool(x,y,n,alpha)作用類似于polyfit.他是一個交互式畫面。eg.y=[1035624108410521015106670496099010508391030985855];x=[6.00002.50007.50008.500010.00007.00003.000011.50005.50006.50004.00009.000011.000012.5000];plot(x,y,'+'),pausex2=x.^2;X=[ones(14,1)x'x2'];[b,bi,r,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