matlab在統(tǒng)計中的應用

1、數(shù)理統(tǒng)計建模Matlab在統(tǒng)計中的應用,山西財經(jīng)大學應用數(shù)學學院 高崇山,一、 概率分布及有關函數(shù),調(diào)用格式為：分布命令符功能命令符(x,參數(shù))[m,v]=分布stat(x,參數(shù))rand產(chǎn)生[0,1]上的隨機數(shù)，randn產(chǎn)生標準正態(tài)分布隨機數(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)在對應于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)生一個m × n的服從參數(shù)為lamda的指數(shù)分布的隨機矩陣,二、 描述性統(tǒng)計 描述性統(tǒng)計就是搜集、整理、加工和分析統(tǒng)計數(shù)據(jù)，使之系統(tǒng)化，以顯示出數(shù)據(jù)的趨勢、特征和數(shù)量關系。,2.1 樣本均值mean和中值median 它們都是樣本數(shù)據(jù)在數(shù)據(jù)分布線上中心位置的度量.,A=[1 2 4 4;3 4 6 6;5 6 8 8;5 6 8 8];mean

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

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

6、的觀察值(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ù)。 百分

7、位數(shù)是用于反映樣本數(shù)據(jù)形態(tài)信息的數(shù)據(jù)統(tǒng)計量，它也可以刻劃數(shù)據(jù)的位置和散布特征。 Y = prctile(X,p) 返回樣本X中大于p%(0<p<100)的值。如果X是向量，則返回X中p百分位數(shù)，若X為矩陣，則返回一個關于每列元素的p百分位數(shù)行向量。注意：p也可以是一個向量，此時，返回一組百分位數(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 相關系數(shù) 相關系數(shù)反映兩個隨即變量之間線性相依程度的變量。R=corrcoef(X)R=corrcoef(x,y)含義同協(xié)方差。同協(xié)方差之間有如下關系,2.5 樣本峰度和偏度偏度描述的是分布的對稱性，其定義為： 當f >0時，表

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

11、i,sigmaci]=normfit(x,alpha)說明：x是樣本(矩陣或向量)；alpha是顯著性水平(默認值為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為

12、可選項，表示顯著性水平；pl僅用于二項分布，表示試驗的次數(shù)；phat為返回的點估計值；pci為返回相應置信區(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),四、 假設經(jīng)驗4.1 單個樣本的t檢驗功能：進行樣本均值的t檢驗格式：h=ttest(x,m)；% 在0.05的顯著性水平下進行t檢驗，以確定在標準差未知的情況下取自正態(tài)分布的均值是否為m，若輸出h=0，則接受零假設，h=1則否定零假設；h=ttest(x,m,alpha)； % alpha 為給定的顯著性水平。[h,sig,ci]=ttest(x,m,alp

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

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

16、test2(x,y,alpha,tail),五、 統(tǒng)計繪圖5.1 box圖boxplot(X,notch,’sym’,vert,whis)X：樣本數(shù)據(jù)；notch=1有切口, notch=0無切口,(默認notch=0)；‘sym’野值標記符號,默認‘+’；vert=0,box圖是水平放置,vert=1是垂直放置(默認)。whis 定義虛線的長度，一般用缺省值.,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)圖見下頁,說明：1.盒子的上下兩條線分別為樣本的75%和25%分位線，中間為樣本中位數(shù)；2.虛線表示樣本的其余部分，位于盒子的上下兩側(cè)；3.‘+’表示野值（奇異值），位于虛線的上方和下方；4.‘切口’表示樣本中位數(shù)的置信區(qū)間。默認狀態(tài)下無切口。,5.2 正態(tài)概率圖正態(tài)概率圖用于判斷

18、樣本數(shù)據(jù)是否服從正態(tài)分布。格式：normplot(X)X：數(shù)據(jù).若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ù)圖用于比較

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

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

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

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

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

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

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

26、值，第三個是F(1,N-2)分布大于F值的概率p，p<alpha時，回歸模型有效。,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(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(x=50)區(qū)間（t分布）[y0 y1 y2]d1=norminv(0.975)*sqrt(s2);y3=y0-d1;y4=y0+d1;[y0 y3 y4] % 預測y(x=50)區(qū)間（N分布）,7.2 多項式回歸1)多項式曲線擬合 多項式回歸的模型為：p(x)=p1

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

31、的是1-alpha的置信區(qū)間y±delta.說明：命令polytool(x,y,n,alpha)作用類似于polyfit.他是一個交互式畫面。,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 多元二項式回歸模型： y = b0 + b1x1+b2x2+b3

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

35、 alpha為可選項。% model=‘interaction’表示回歸模型包括常數(shù)項、線性項、交叉項；% model=‘quadratic’表示回歸模型包括常數(shù)項、線性項、交叉項、二次項；% model=‘purequadratic’表示回歸模型包括常數(shù)項、線性項、二次項；rstool(x,y,’model’,alpha,’xname’,’yname’)% 可以在x軸和y軸上分別標上相應的名稱。,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);