數理統計建模

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

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

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

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

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

輸入y(因變量，列向量)；x(第一列全為1，第二列為x1的觀察值，第三列為x2的觀察值，…)；alpha為顯著性水平α(默認值為0.05)；輸出b為(β0,β2,…β2,βm)的估計值；bint是(β0,β2,…β2,βm)的置信區間；r是殘差(觀察值與預測值之差，為列向量，主要用于探測模型假設的合理性)，rint是殘差的置信區間；s包含3個統計量：第一個是決定系數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(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(x=50)區間〔t分布〕[y0y1y2]d1=norminv(0.975)*sqrt(s2);y3=y0-d1;y4=y0+d1;[y0y3y4]%預測y(x=50)區間〔N分布〕7.2多項式回歸1)多項式曲線擬合多項式回歸的模型為：p(x)=p1xn+p2xn-2+…+pnx+pn+1格式為：[p,s]=polyfit(x,y,n)其中：n為擬合次數；x,y分別是自變量和因變量；s是一個矩陣，用于polyval函數，可進行預測的誤差估計；p為系數向量(p1,p2,…,pn,pn+1)的估計值。2〕多項式預測和置信區間的評估[y,delta]=polyconf(p,x,s,alpha)其中：p,s是擬合輸出的結果，x是要預測的點，alpha是置信度，輸出的是1-alpha的置信區間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