第10章 用MATLAB解決概率問題

MATLAB是一種功能強大、運算效率極高的數(shù)值計算軟件。最初，它是一種專門用于矩陣運算的軟件，經(jīng)過多年發(fā)展，MATLAB已經(jīng)成為一種功能強大的軟件，幾乎可以解決科學計算中的任何問題。本章將利用Matlab來解決概率統(tǒng)計學中的概率分布、數(shù)字特征、參數(shù)估計以及假設(shè)檢驗等問題。第10章用Matlab解決概率問題數(shù)據(jù)分析離散型隨機變量的概率及概率分布連續(xù)型隨機變量的概率及其分布數(shù)字特征二維隨機向量的數(shù)字特征統(tǒng)計直方圖參數(shù)估計假設(shè)檢驗方差分析與回歸分析10.2、離散型隨機變量的概率及概率分布(1)分布律二項分布的概率值

格式

binopdf(k,n,p)

說明

n:試驗總次數(shù)；p:每次試驗事件A發(fā)生的概率；k:事件A發(fā)生k次。泊松分布的概率值

格式

poisspdf(k,lambda)

說明

k:事件A發(fā)生k次;lambda:參數(shù)超幾何分布的概率值

格式

hygpdf(K,N,M,n)

說明

K:抽得次品數(shù)；N：產(chǎn)品總數(shù)；M：次品總數(shù)；n:抽取總數(shù).(2)累積概率值(隨機變量X<K的概率之和)

二項分布的累積概率值

格式binocdf(k,n,p)

說明

n:試驗總次數(shù)；p:每次試驗事件A發(fā)生的概率；k:事件A發(fā)生k次。泊松分布的累積概率值

格式poisscdf(k,lambda)

說明k:事件A發(fā)生k次;lambda:參數(shù)超幾何分布的累積概率值格式hygcdf(K,N,M,n)

說明

K:抽得次品數(shù)；N：產(chǎn)品總數(shù)；M：次品總數(shù)；n:抽取總數(shù).應用舉例例1

某機床出次品的概率為0.01,求生產(chǎn)1000件產(chǎn)品中：（1）恰有一件次品的概率；（2）至少有一件次品的概率。

解：此問可看作是1000次獨立重復試驗，每次試驗出次品的概率為0.01，恰有一件次品的概率,在Matlab命令窗口鍵入：

>>p=binopdf(1,1000,0.01)

顯示結(jié)果為：

p=0.3681（2）至少有一件次品的概率,在Matlab命令窗口鍵入：

>>p=1-binocdf(1,1000,0.01)

顯示結(jié)果為：p=0.6323例2自1875年到1955年中的某63年間，某城市夏季（5-9月間）共發(fā)生暴雨160次，試求在一個夏季中發(fā)生k次（k=0,1,2,…,8）暴雨的概率（設(shè)每次暴雨以1天計算）。解：一年夏天共有天數(shù)為

n=31+30+31+31+30=153

故可知夏天每天發(fā)生暴雨的概率約為很小，n=153較大，可用泊松分布近似。應用舉例10.3連續(xù)型隨機變量的概率及其分布(1)概率密度函數(shù)值利用專用函數(shù)計算概率密度函數(shù)值，如下表。分布調(diào)用函數(shù)均勻分布unifpdf(x,a,b)指數(shù)分布exppdf(x,lambda)正態(tài)分布normpdf(x,mu,sigma)分布chi2pdf(x,n)T分布tpdf(x,n)F分布fpdf(x,n1,n2)應用舉例例5計算正態(tài)分布N（0，1）下的在點0.6633的值。在Matlab命令窗口鍵入：

>>normpdf(0.6633,0,1)回車后顯示結(jié)果為：ans=0.3202舉例應用例6繪制卡方分布密度函數(shù)在n分別等于5，5，20時的圖形程序：x=0:0.1:30;y1=chi2pdf(x,5);plot(x,y1,':')holdon%保留當前圖形y2=chi2pdf(x,15);plot(x,y2,'+')y3=chi2pdf(x,20);plot(x,y3,'o')axis([0,30,0,0.2])%控制圖形在坐標軸上的范圍xlabel(‘圖2-1’)%給軸標注“圖2-1”結(jié)果為下圖(2)分布函數(shù)利用專用函數(shù)計算累積概率函數(shù)值，即

常用專用函數(shù)如下表。分布調(diào)用函數(shù)均勻分布unifcdf(x,a,b)指數(shù)分布expcdf(x,lambda)正態(tài)分布normcdf(x,mu,sigma)卡方分布chi2cdf(x,n)T分布tcdf(x,n)F分布fcdf(x,n1,n2)應用舉例例7某公共汽車站從上午7：00起每15分鐘來一班車。若某乘客在7：00到7：30間任何時刻到達此站是等可能的，試求他候車的時間不到5分鐘的概率。解：設(shè)乘客7點過X分鐘到達此站，則X在[0，30]內(nèi)服從均勻分布，當且僅當他在時間間隔（7：10，7：15）或（7：25，7：30）內(nèi)到達車站時，候車時間不到5分鐘。故其概率為：P1=P{10<X<15}+P{25<X<30}程序：

>>

formatrat>>p1=unifcdf(15,0,30)-unifcdf(10,0,30);>>p2=unifcdf(30,0,30)-unifcdf(25,0,30);>>p=p1+p2則結(jié)果顯示為：p=1/3應用舉例例9設(shè)隨機變量X的概率密度為

確定常數(shù)c;求X落在區(qū)間（-1/2，1/2）內(nèi)的概率；求X的分布函數(shù)F（x）程序（1）：>>symscx>>px=c/sqrt(1-x.^2);>>Fx=int(px,x,-1,1)則結(jié)果顯示如下：Fx=pi*c由pi*c=1得

c=1/pi程序（2）：>>symsx>>c='1/pi';>>px=c/sqrt(1-x.^2);>>format>>p1=int(px,x,-1/2,1/2)

則結(jié)果顯示如下：

p1=1/3程序(3)>>symsxt>>c='1/pi';>>px=c/sqrt(1-t.^2);>>format>>Fx=int(px,t,-1,x)

則結(jié)果顯示如下：Fx=1/2*(2*asin(x)+pi)/pi

所以X的分布函數(shù)為：逆累積概率值已知，求x。x為臨界值，常用臨界值如表調(diào)用函數(shù)注釋X=unifinv(p,a,b)[a,b]上均勻分布逆累積分布函數(shù)X=expinv(p,mu)指數(shù)逆累積分布函數(shù)X=norminv(p,mu,sigma)正態(tài)逆累積分布函數(shù)X=chi2inv(p,n)卡方逆累積分布函數(shù)X=tinv(p,n)T分布逆累積分布函數(shù)X=finv(p,n1,n2)F分布逆累積分布函數(shù)應用舉例例13公共汽車門的高度是按成年男子與車門頂碰頭的機會不超過2%設(shè)計的。設(shè)男子身高X（單位：cm）服從正態(tài)分布N（168，7），求車門的最低高度。解：設(shè)h為車門高度，X為身高，求滿足條件由已知，P{X>=x}=0.02,即P{X<=x}=0.98>>norminv(0.98,168,7)ans=182.3762所以至少為182.4厘米應用舉例例14設(shè)二維隨機變量（X，Y）的聯(lián)合密度為：

求（1）P{0<X<1,0<Y<1};(2)(X,Y)落在x+y=1,x=0,y=0所圍成的區(qū)域內(nèi)的概率。程序：>>symsxy>>f=exp(-x-y);>>P_XY=int(int(f,y,0,1),x,0,1)>>P_G=int(int(f,y,0,1-x),x,0,1)運行結(jié)果顯示如下：P_XY=exp(-2)-2*exp(-1)+1P_G=-2*exp(-1)+110.4數(shù)字特征(1)數(shù)學期望離散型隨機變量X的期望計算求和函數(shù)：sum(X)

說明：若X為向量，則sum（X）為X中的各元素之和，返回一個數(shù)值；若X為矩陣，則sum（X）為X中各列元素之和，返回一個行向量。求均值函數(shù)：mean(X)

說明：若X為向量，則sum（X）為X中的各元素的算術(shù)平均值，返回一個數(shù)值；若X為矩陣，則sum（X）為X中各列元素的算術(shù)平均值，返回一個行向量例16隨機抽取6個滾珠測得直徑（mm）如下：

11.7012.2111.9011.9112.3212.32

試求樣本平均值。程序：>>X=[11.7012.2111.9011.9112.3212.32];>>mean(X)則結(jié)果顯示如下：

ans=12.0600

應用舉例連續(xù)型隨機變量的期望應用舉例例17已知隨機變量X的概率

求EX和E(4X-1)。程序：解：在Matlab編輯器中建立M文件LX0817.m：symsxp_x=3*x^2;EX=int(x*p_x,0,1)EY=int((4*x-1)*p_x,0,1)運行結(jié)果為：EX=3/4EY=2(2)方差離散型隨機變量的方差及樣本方差方差設(shè)X的分布律為由則方差DX=sum(X.^2*P)-(EX).^2標準差：應用舉例例19設(shè)隨機變量X的分布律為：X-2-1012P0．30．10．20．10．3求D（X），D（X^2-1）。程序：>>X=[-2-1012];>>p=[0.30.10.20.10.3];>>EX=sum(X.*p)>>Y=X.^2-1>>EY=sum(Y.*p)>>DX=sum(X.^2.*p)-EX.^2>>DY=sum(Y.^2.*p)-EY.^2運行后結(jié)果顯示如下：EX=0Y=30-103EY=1.6000DX=2.6000DY=3.0400連續(xù)型隨機變量的方差利用求解。例21設(shè)X的概率密度為：求DX，D（2X+1）解：程序：>>symsx>>px=1./(pi*sqrt(1-x.^2));>>EX=int(x*px,-1,1)>>Dx=int(x.^2.*px,-1,1)>>y=2*x+1;>>EY=int(y.*px,-1,1)>>DY=int(y.^2.*px,-1,1)-EY.^2運行結(jié)果顯示如下：EX=0DX=1/2EY=1DY=2（3）常用分布的期望與方差求法在統(tǒng)計工具箱中，用’stat’結(jié)尾的函數(shù)可以計算給定參數(shù)的某種分布的均值和方差。調(diào)用函數(shù)參數(shù)說明注釋[M，V]=binostat(N,p)N：試驗總次數(shù)P:二項分布的概率二項分布的期望和方差[M,V]=poisstat(Lambda)Lambda:泊松分布的參數(shù)泊松分布的期望和方差[M，V]=hygestat(M,K,N)M,K,N:超幾何分布的參數(shù)超幾何分布期望和方差[M,V]=unifstat(a,b)a,b:均勻分布區(qū)間端點均勻分布的期望和方差[M,V]=normstat(mu,sigma)Mu:正態(tài)分布的均值Sigma：標準差正態(tài)分布的期望和方差調(diào)用函數(shù)參數(shù)說明注釋[M,V]=tstat(nu)nu:T分布的參數(shù)T分布的期望和方差[M,V]=chi2stat(nu)nu:卡方分布的參數(shù)卡方分布的期望和方差[M,V]=fstat(n1,n2)n1,n2:F分布2個自由度F分布的期望和方差[M,V]=geostat(p)p:幾何分布的概率參數(shù)幾何分布的期望和方差[M,V]=expstat(mu)mu:指數(shù)分布的參數(shù)指數(shù)分布的期望和方差應用舉例例24求參數(shù)為6的泊松分布的期望和方差。程序：>>[M,V]=poisstat(6)則結(jié)果顯示如下：

M=6V=610.5二維隨機變量的數(shù)字特征(1)期望

根據(jù)二維隨機變量期望的定義構(gòu)造函數(shù)計算。下面分別就離散和連續(xù)的情況舉例說明。應用舉例例4.1設(shè)（X，Y）的聯(lián)合分布為Z=X-Y，求EZ。XY-112-15/202/206/2023/203/201/20程序：>>X=[-12];>>Y=[-112];>>fori=1:2forj=1:3Z(i,j)=X(i)-Y(j);endend>>P=[5/202/206/20;3/203/201/20];>>EZ=sum(sum(Z.*P))%將Z與P對應相乘相加運行結(jié)果顯示如下：

EZ=-0.5000應用舉例例26射擊試驗中，在靶平面建立以靶心為原點的直角坐標系，設(shè)X，Y分別為彈著點的橫坐標和縱坐標，它們相互獨立且均服從N(0,1)，求彈著點到靶心距離的均值。解：設(shè)彈著點到靶心距離為，則求EZ。聯(lián)合概率密度：期望為：程序：>>symsxyrt>>pxy=1/(2*pi)*exp(-1/2*(x.^2+y.^2));>>EZ=int(int(r*1/(2*pi)*exp(-1/2*r.^2)*r,r,0,+inf),t,0,2*pi)運行結(jié)果：EZ=

1/2*2^(1/2)*pi^(1/2)即（2）協(xié)方差Matlab提供了求協(xié)方差的函數(shù)：cov(X)

％X為向量時，返回此向量的方差；X為矩陣時，返回此矩陣的協(xié)方差矩陣cov(X，Y)

%返回X與Y的協(xié)方差，X與Y同維數(shù)cov(X，0)

%返回X的樣本協(xié)方差，置前因子為1/（n-1）cov(X，1)

%返回X的協(xié)方差，置前因子為1/n應用舉例例27設(shè)(X,Y)的聯(lián)合概率密度為：求DX，DY和解：程序：>>symsxy>>pxy=1/8*(x+y);>>EX=int(int(x*pxy,y,0,2),0,2)>>EY=int(int(y*pxy,x,0,2),0,2)>>EXX=int(int(x^2*pxy,y,0,2),0,2)>>EYY=int(int(y^2*pxy,x,0,2),0,2)>>EXY=int(int(x*y*pxy,x,0,2),0,2)>>DX=EXX-EX^2>>DY=EYY-EY^2>>DXY=EXY-EX*EY運行結(jié)果顯示如下：EX=7/6EY=7/6EXX=5/3EYY=5/3EXY=4/3DX=11/36DY=11/36DXY=-1/36例29求一個隨機矩陣的協(xié)方差。在命令窗口鍵入>>d=rand(2,4)則結(jié)果為d=

0.95010.60680.89130.45650.23110.48600.76210.0185>>cov1=cov(d)cov1=0.25850.04340.04640.15740.04340.00730.00780.02650.04640.00780.00830.02830.15740.02650.02830.0959（3）相關(guān)系數(shù)Matlab提供了求相關(guān)系數(shù)的函數(shù)。corrcoef(X,Y)%返回列向量X，Y的相關(guān)系數(shù)corrcoef(X)%返回矩陣X的列向量的相關(guān)系數(shù)矩陣例4.6設(shè)（X，Y）在單位圓上服從均勻分布，即有聯(lián)合密度求，，及。解：程序：>>symsxyrt>>PXY='1/pi';>>EX=int(int(r^2*cos(t)*PXY,r,0,1),0,2*pi)>>EY=int(int(r^2*sin(t)*PXY,r,0,1),0,2*pi)>>EXX=int(int(r^3*cos(t)^2*PXY,r,0,1),0,2*pi)>>EYY=int(int(r^3*sin(t)^2*PXY,r,0,1),0,2*pi)>>EXY=int(int(r^3*sin(t)*cos(t)*PXY,r,0,1),0,2*pi)>>DXX=EXX-(EX)^2>>DYY=EYY-(EY)^2>>DXY=EXY-EX*EY>>ro_XY=DXY/sqrt(DXX*DYY)運行后結(jié)果顯示如下：EX=0EY=0EXX=1/4EYY=1/4EXY=0DXX=1/4DYY=1/4DXY=0Ro_XY=010.6統(tǒng)計直方圖函數(shù)

hist（Z，n）

%直角坐標系下的統(tǒng)計直方圖

n表示直方圖的區(qū)間數(shù)，缺省時n=10函數(shù)rose(theta,n)

%極坐標系下角度直方圖

n是在[0，2]范圍內(nèi)所分區(qū)域數(shù)，缺省時n=20,theta為指定的弧度數(shù)據(jù)。應用舉例例32某食品廠為加強質(zhì)量管理，對生產(chǎn)的罐頭重量X進行測試，在某天生產(chǎn)的罐頭中抽取了100個，其重量測試數(shù)據(jù)記錄如下：

332330338336333332336331334338346346340344342344345340344344

343344342343345339350337345349336348344345332342342340350343347340344353340340356346345346340339342352342350348344350335340338345345349336342338343343

341347341347344339347348343347346344345350341338343339343346342339343350341346341345344342試根據(jù)以上數(shù)據(jù)作出X的頻率直方圖。程序：>>在Matlab編輯器中建立M文件LX0832.m：X=[332330338336333332336331334338...346346340344342344345340344344...343344342343345339350337345349...336348344345332342342340350343...347340344353340340356346345346...340339342352342350348344350335...340338345345349336342338343343...341347341347344339347348343347...346344345350341338343339343346...342339343350341346341345344342];>>hist(X,13)該例的統(tǒng)計直方圖如下：10.7參數(shù)估計樣本數(shù)字特征法樣本均值：

mx=1/n*sum(x)

樣本方差：

ss=1/(n-1)*sum((x-mx).^2)（1）點估計（2）最大似然估計調(diào)用函數(shù)函數(shù)說明binofit(X,N)[PHAT,PCI]=binofit(X,N,ALPHA)水平的參數(shù)估計和置信區(qū)間poissfit(X)[LAMBDHAT,LAMBDPCI]=poissfit(X,ALPHA)參數(shù)估計和置信區(qū)間normfit(X,ALPHA)[MUHAT,SIGMAHAT,MUPCI,SIGMACI]=normfit(X,ALPHA)水平的期望、方差值和置信區(qū)間unfit(X,ALPHA)[AHAT,BHAT,ACI,BCI]=unfit(X,ALPHA)水平的參數(shù)估計和置信區(qū)間expfit(X)[MUHHAT,MUCI]=expfit(X,ALPHA)水平的參數(shù)估計和置信區(qū)間應用舉例例34設(shè)某種清漆的9個樣品，其干燥時間（以小時計）分別為

5.04.74.85.56.05.34.65.14.0

設(shè)干燥時間總體服從正態(tài)分布N（）求和的置信度為0.95的置信區(qū)間（未知）程序：>>X=[5.04.74.85.56.05.34.65.14.0

];>>[MUHAT,SIGMAHAT,MUPCI,SIGMACI]=normfit(X,0.05)運行后結(jié)果顯示如下：

muhat=5%的最大似然估計值sigmahat=0.5745%的最大似然估計值muci=4.55845.4416%的置信區(qū)間sigmaci=0.38801.1005%的置信區(qū)間10.8假設(shè)檢驗在Matlab中，假設(shè)檢驗問題都提出兩種假設(shè)：即原假設(shè)和備擇假設(shè)。對于正態(tài)總體均值的假設(shè)檢驗給出了檢驗函數(shù)：ztest

已知，檢驗正態(tài)總體均值；ttest

未知，檢驗正態(tài)總體均值；ttest2

兩個正態(tài)總體均值比較。對于一般連續(xù)型總體一致性的檢驗，給出了檢驗方法—秩和檢驗，由函數(shù)ranksum實現(xiàn)。1單個正態(tài)總體N（）的假設(shè)檢驗

已知，對期望的假設(shè)檢驗—Z檢驗法調(diào)用函數(shù)H=ztest（X，m,sigma）

H=ztest（X，m,sigma,alpha）

[H,sig,ci]=ztest（X，m,sigma,alpha,tail）說明：X：樣本；m:期望值；sigma:正態(tài)總體標準差；alpha:經(jīng)驗水平；

tail:備擇假設(shè)的選項，若tail=0(缺省)，則；若tail=1，則；若tail=-1，則。即tail=0(缺省)為雙邊檢驗，其余為單邊檢驗問題。

H：檢驗結(jié)果，分兩種情況：若H=0，則在水平下，接受原假設(shè)；若H=1，則在水平下，拒絕原假設(shè)。

sig為當原假設(shè)為真時（即成立），得到觀察值的概率，當sig為小概率時，則對原假設(shè)提出質(zhì)疑。Ci:均值的1-alpha置信區(qū)間。應用舉例例36某車間用一臺包裝機包裝葡萄糖，包得的袋裝糖重是一個隨機變量，它服從正態(tài)分布。當機器正常時，其均值為0.5公斤，標準差為0.015。某日開工后檢驗包裝機是否正常，隨機地抽取所包裝的糖9袋，稱得凈重為：（公斤）

0.4970.5060.5180.5240.4980.5110.520.5150.512問機器是否正常？解：已知，在水平=0.05下檢驗假設(shè)：原假設(shè)：備擇假設(shè)：

程序：>>X=[0.4970.5060.5180.5240.4980.5110.520.5150.512

];>>[H,SIG]=ztest(X,0.5,0.015,0.05,0)運行后顯示結(jié)果如下：

H=1SIG=0.0248結(jié)果表明：H=1，說明在水平=0.05下，可拒絕原假設(shè)，即認為包裝機工作不正常。1單個正態(tài)總體N（）的假設(shè)檢驗

未知，對期望的假設(shè)檢驗—t檢驗法調(diào)用函數(shù)H=ttest(X,m,sigma)%在水平=sigma下檢驗是否成立。說明：X：樣本；m:期望值；alpha:經(jīng)驗水平；

tail:備擇假設(shè)的選項，若tail=0(缺省)，則備擇假設(shè)為

；若tail=1，則；若tail=-1，則。即tail=0(缺省)為雙邊檢驗，其余為單邊檢驗問題。

H：檢驗結(jié)果，分兩種情況：若H=0，則在水平下，接受原假設(shè)；若H=1，則在水平下，拒絕原假設(shè)。

sig為當原假設(shè)為真時（即成立），得到觀察值的概率，當sig為小概率時，則對原假設(shè)提出質(zhì)疑。Ci:均值的1-alpha置信區(qū)間。應用舉例例37某種電子元件的壽命X（以小時計）服從正態(tài)分布，均未知，現(xiàn)測得16只元件壽命如下：15928010121222437917