概率統(tǒng)計(jì)的數(shù)值實(shí)驗(yàn)-MATLAB在概率統(tǒng)計(jì)教學(xué)中的應(yīng)用

引言

而MATLAB軟件具有簡(jiǎn)單易學(xué)、易操作和繪圖功能強(qiáng)等特點(diǎn)，利用MATLAB軟件的圖形可視功能將概率統(tǒng)計(jì)的內(nèi)容用圖形表示出來(lái)，通過圖形讓學(xué)生加深理解，以達(dá)到事半功倍的效果。

概率論與數(shù)理統(tǒng)計(jì)知識(shí)比較抽象，邏輯性較強(qiáng)。因此，建議讓學(xué)生結(jié)合理論和公式推導(dǎo)，進(jìn)行數(shù)值試驗(yàn)和相關(guān)調(diào)查，直觀地感受數(shù)學(xué)概念和理論，從而提高學(xué)生解決實(shí)際問題的信心和能力。概率論

1.rand(m,n)：生成m×n的隨機(jī)矩陣，每個(gè)元素都在(0,1)

間，生成方式為均勻分布。

2.randn(m,n)：生成m×n的隨機(jī)矩陣，每個(gè)元素都在(0,1)

間，生成方式為正態(tài)分布。

3.randperm(m)：生成一個(gè)1~m的隨機(jī)整數(shù)排列。

4.perms(1:n)：生成一個(gè)1~n的全排列，共n!個(gè)。

5.取整函數(shù)系列：

（1）fix(x)：截尾法取整；

（2）floor(x)：退一法取整（不超過x的最大整數(shù)）；

（3）ceil(x)：進(jìn)一法取整（=floor(x)+1）；

（4）round(x)：四舍五入法取整。

6.unique(a)：合并a中相同的項(xiàng)。

7.prod(x)：向量x的所有分量元素的積。一、MATLAB常用的與隨機(jī)數(shù)產(chǎn)生相關(guān)的函數(shù)：示例：

>>rand(1)%生成一個(gè)(0,1)間的隨機(jī)數(shù)

ans=0.8147

>>rand(2,2)%生成一個(gè)2×2階(0,1)間的隨機(jī)數(shù)矩陣

ans=0.91340.0975

0.63240.2785

>>randperm(5)%生成一個(gè)1~5的隨機(jī)整數(shù)排列

ans=41523

>>a=[1242332];

unique(a)

ans=1234

例1

隨機(jī)投擲均勻硬幣，觀察國(guó)徽朝上與國(guó)徽

朝下的頻率。解>>n=3000~100000000;m=0;fori=1:nt=randperm(2);%生成一個(gè)1~2的隨機(jī)整數(shù)排列x=t-1;%生成一個(gè)0~1的隨機(jī)整數(shù)排列y=x(1);ify==0;m=m+1;endendp1=m/np2=1-p1

試驗(yàn)次數(shù)n300050001萬(wàn)2萬(wàn)3萬(wàn)國(guó)徽朝上頻率0.50400.50060.48790.49990.5046國(guó)徽朝下頻率0.49600.49940.51210.50010.4954試驗(yàn)次數(shù)n5萬(wàn)10萬(wàn)100萬(wàn)100萬(wàn)1億國(guó)徽朝上頻率0.50210.49990.49990.50010.5000國(guó)徽朝下頻率0.49790.50010.50010.49990.5000可見當(dāng)時(shí)，解記事件為第i個(gè)人拿到自已槍，事件為第i個(gè)人沒拿到自己槍，易知：又記為沒有一個(gè)人拿到自己槍的概率。

有乘法公式可知：例2

某班有n個(gè)人，每人各有一支槍，這些槍外形一樣。某次夜間緊急集合，若每人隨機(jī)地取走一支槍，問沒有一個(gè)人拿到自己槍的概率是多少？于是

所以

特別地，當(dāng)n較大時(shí)，。因此，可隨機(jī)模擬出沒有人拿到自己槍的頻率，根據(jù)頻率的穩(wěn)定性，近似當(dāng)做概率，然后去估計(jì)自然常數(shù)e。算法如下：

1、產(chǎn)生n個(gè)隨機(jī)數(shù)的隨機(jī)序列；

2、檢驗(yàn)隨機(jī)列與自然列是否至少有一個(gè)配對(duì)；

3、對(duì)沒有一個(gè)配對(duì)的序列進(jìn)行累積p；

4、重復(fù)1、2、3步m

次；

5、估計(jì)。

具體程序及相關(guān)結(jié)果為（注：自然常數(shù)

e≈2.7183）：>>m=40000;n=50;p=0;forj=1:mk=0;sui=randperm(n);fori=1:nifsui(i)==ik=k+1;elsek=k;endendifk==0p=p+1;elsep=p;endende=m/pe=2.7313模擬次數(shù)m400004000040000人數(shù)n100020005000e2.71552.70822.7202模擬次數(shù)m400040000400000人數(shù)n505050e2.73792.73132.7194

設(shè)針與平行線的夾角為，針的中心與最近直線的距離為。針與平行線相交的充要條件是，則所求概率為

故可得的近似計(jì)算公式，其中n為隨機(jī)試驗(yàn)次數(shù)，m為針與平行線相交的次數(shù)。例3Buffon投針實(shí)驗(yàn)

在畫有許多間距為的等距平行線的白紙上，隨機(jī)投擲一根長(zhǎng)為的均勻直針，求針與平行線相交的概率，并計(jì)算的近似值。解

>>clear,clfn=10000000;l=0.5;m=0;d=1;fori=1:nx=l/2*sin(rand(1)*pi);y=rand(1)*d/2;ifx>=ym=m+1;endendp1=m/npai=2*n*l/(m*d)試驗(yàn)次數(shù)n5千1萬(wàn)10萬(wàn)100萬(wàn)1000萬(wàn)針長(zhǎng)l/平行間距d3/103/103/103/103/10相交頻率0.18360.19710.18870.19050.1912π的近似值3.26803.04413.17983.14983.1387試驗(yàn)次數(shù)n5千1萬(wàn)10萬(wàn)100萬(wàn)1000萬(wàn)針長(zhǎng)l/平行間距d2/52/52/52/52/5相交頻率0.24960.25620.25490.25440.2543π的近似值3.20513.12263.13863.14513.1433試驗(yàn)次數(shù)n5千1萬(wàn)10萬(wàn)100萬(wàn)1000萬(wàn)針長(zhǎng)l/平行間距d1/21/21/21/21/2相交頻率0.32540.31480.31580.31780.3183π的近似值3.07313.17663.16673.14703.1417試驗(yàn)次數(shù)n5千1萬(wàn)10萬(wàn)100萬(wàn)1000萬(wàn)針長(zhǎng)l/平行間距d4/54/54/54/54/5相交頻率0.51420.51340.50860.50930.5093π的近似值3.11163.11653.14603.14183.1418試驗(yàn)次數(shù)n5千1萬(wàn)10萬(wàn)100萬(wàn)1000萬(wàn)針長(zhǎng)l/平行間距d17/2017/2017/2017/2017/20相交頻率0.54320.54520.54200.54120.5410π的近似值3.12963.11813.13663.14133.1426試驗(yàn)次數(shù)n5千1萬(wàn)10萬(wàn)100萬(wàn)1000萬(wàn)針長(zhǎng)l/平行間距d9/109/109/109/109/10相交頻率0.58600.57000.57560.57330.5731π的近似值3.07173.15793.12723.13953.1410例4在100個(gè)人的團(tuán)體中，不考慮年齡差異，研究是否有兩個(gè)以上的人生日相同。假設(shè)每人的生日在一年365天中的任意一天是等可能的，那么隨機(jī)找n個(gè)人(不超過365人)。

(1)求這n個(gè)人生日各不相同的概率是多少？從而求這n個(gè)人中至少有兩個(gè)人生日相同這一隨機(jī)事件發(fā)生的概率是多少？

(2)近似計(jì)算在30名學(xué)生的一個(gè)班中至少有兩個(gè)人生日相同的概率是多少？解:(1)>>clear,clfforn=1:100p0(n)=prod(365:-1:365-n+1)/365^n;p1(n)=1-p0(n);endp1=ones(1,100)-p0;n=1:100;plot(n,p0,n,p1,'--')xlabel('人數(shù)'),ylabel('概率')legend('生日各不相同的概率','至少兩人生日相同的概率')axis([0100-0.11.199]),gridonp1(30)=0.7063,p1(60)=0.9941

分析：在30名學(xué)生中至少兩人生日相同的概率為70.63％。下面進(jìn)行計(jì)算機(jī)仿真。

隨機(jī)產(chǎn)生30個(gè)正整數(shù)，代表一個(gè)班30名學(xué)生的生日，然后觀察是否有兩人以上生日相同。當(dāng)30個(gè)人中有兩人生日相同時(shí)，輸出“1”，否則輸出“0”。如此重復(fù)觀察100次，計(jì)算出這一事件發(fā)生的頻率。

(2)>>clear,clfn=0;form=1:100%做100次隨機(jī)試驗(yàn)

y=0;x=1+fix(365*rand(1,30));%產(chǎn)生30個(gè)隨機(jī)數(shù)

fori=1:29%用二重循環(huán)尋找30個(gè)隨機(jī)數(shù)中是否有相同數(shù)

forj=i+1:30ifx(i)==x(j)y=1;break;endendendn=n+y;%累計(jì)有兩人生日相同的試驗(yàn)次數(shù)endf=n/m%計(jì)算頻率f=0.6900f=0.7900f=0.6700f=0.7300f=0.7500f=0.6900f=0.7200f=0.6700f=0.6800……重復(fù)觀察，數(shù)據(jù)如下：例5Galton釘板模型和二項(xiàng)分布

Galton釘板試驗(yàn)是由英國(guó)生物統(tǒng)計(jì)學(xué)家和人類學(xué)家Galton設(shè)計(jì)的。故而得名。

通過模擬Calton釘板試驗(yàn)，觀察和體會(huì)二項(xiàng)分布概率分布列的意義、形象地理解DeMoivre-Laplace中心極限定理。共15層小釘Ox-8-7-6-5-4-3-2-1

12345678高爾頓釘板試驗(yàn)小球最后落入的格數(shù)?記小球向右落下的次數(shù)為

則記小球向左落下的次數(shù)為

則符號(hào)函數(shù),大于0返回1,小于0返回-1,等于0返回0

高爾頓(FrancisGalton,1822-1911)英國(guó)人類學(xué)家和氣象學(xué)家Ox-8-7

-6

-5

-4

-3

-2

-1

12345678記則近似高爾頓釘板試驗(yàn)？什么曲線共15層小釘小球碰第層釘后向右落下小球碰第層釘后向左落下

模擬Galton釘板試驗(yàn)的步驟：

(1)確定釘子的位置：將釘子的橫、縱坐標(biāo)存儲(chǔ)在兩個(gè)矩陣X和Y中。

(2)在Galton釘板試驗(yàn)中，小球每碰到釘子下落時(shí)都具有兩種可能性，設(shè)向右的概率為p，向左的概率為q＝1-p，這里p=0.5，表示向左向右的機(jī)會(huì)是相同的。模擬過程如下：首先產(chǎn)生一均勻隨機(jī)數(shù)u，這只需調(diào)用隨機(jī)數(shù)發(fā)生器指令rand(m,n)。

rand(m,n)指令：用來(lái)產(chǎn)生m×n個(gè)(0,1)區(qū)間中的隨機(jī)數(shù)，并將這些隨機(jī)數(shù)存于一個(gè)m×n矩陣中，每次調(diào)用rand(m,n)的結(jié)果都會(huì)不同。如果想保持結(jié)果一致，可與rand(‘seed’,s)配合使用，這里s是一個(gè)正整數(shù)，例如>>rand('seed',1),u=rand(1,6)u=0.51290.46050.35040.09500.43370.7092而且再次運(yùn)行該指令時(shí)結(jié)果保持不變。除非重設(shè)種子seed的值，如>>rand('seed',2),u=rand(1,6)u=0.02580.92100.70080.19010.86730.4185這樣結(jié)果才會(huì)產(chǎn)生變化。

將[0,1]區(qū)間分成兩段，區(qū)間[0,p)和[p,1]。如果隨機(jī)數(shù)u屬于[0,p)，讓小球向右落下；若u屬于[p,1]，讓小球向左落下。將這一過程重復(fù)n次，并用直線連接小球落下時(shí)所經(jīng)過的點(diǎn)，這樣就模擬了小球從頂端隨機(jī)地落人某一格子的過程。

(3)模擬小球堆積的形狀。輸入扔球次數(shù)m(例如m＝50、100、500等等)，計(jì)算落在第i個(gè)格子的小球數(shù)在總球數(shù)m中所占的比例，這樣當(dāng)模擬結(jié)束時(shí)，就得到了頻率

用頻率反映小球的堆積形狀。

(4)用如下動(dòng)畫指令制作動(dòng)畫：

movien(n)：創(chuàng)建動(dòng)畫矩陣；制作動(dòng)畫矩陣數(shù)據(jù)；

Getframe：拷貝動(dòng)畫矩陣；

movie(Mat,m)：播放動(dòng)畫矩陣m次。

M文件如下：解：>>clear,clf,m=100;n=5;y0=2;%設(shè)置參數(shù)ballnum=zeros(1,n+1);p=0.5;q=1-p;fori=n+1:-1:1%創(chuàng)建釘子的坐標(biāo)x，yx(i,1)=0.5*(n-i+1);y(i,1)=(n-i+1)+y0;forj=2:ix(i,j)=x(i,1)+(j-1)*1;y(i,j)=y(i,1);endendmm=moviein(m);%動(dòng)畫開始，模擬小球下落路徑fori=1:ms=rand(1,n);%產(chǎn)生n個(gè)隨機(jī)數(shù)

xi=x(1,1);yi=y(1,1);k=1;l=1;%小球遇到第一個(gè)釘子

forj=1:nplot(x(1:n,:),y(1:n,:),‘o’,x(n+1,:),y(n+1,:),‘.-’),%畫釘子的位置axis([-2n+20y0+n+1]),holdon

k=k+1;%小球下落一格

ifs(j)>pl=l+0;%小球左移

elsel=l+1;%小球右移

endxt=x(k,l);yt=y(k,l);%小球下落點(diǎn)的坐標(biāo)

h=plot([xi,xt],[yi,yt]);axis([-2n+20y0+n+1])%畫小球運(yùn)動(dòng)軌跡

xi=xt;yi=yt;endballnum(l)=ballnum(l)+1;%計(jì)數(shù)

ballnum1=3*ballnum./m;bar([0:n],ballnum1),axis([-2n+20y0+n+1])%畫各格子的頻率

mm(i)=getframe;%存儲(chǔ)動(dòng)畫數(shù)據(jù)

holdoffendmovie(mm,1)%播放動(dòng)畫一次……①概率密度函數(shù)（pdf），求隨機(jī)變量X在x點(diǎn)處的概率密度值②累積分布函數(shù)（cdf），求隨機(jī)變量X在x點(diǎn)處的分布函數(shù)值③逆累積分布函數(shù)（inv），求隨機(jī)變量X在概率點(diǎn)處的分布函數(shù)反函數(shù)值④均值與方差計(jì)算函數(shù)（stat），求給定分布的隨機(jī)變量X的數(shù)學(xué)期望E（X）和方差var（X）。⑤隨機(jī)數(shù)生成函數(shù)（rnd），模擬生成指定分布的樣本數(shù)據(jù)。二、MATLAB為常見自然概率分布提供了下列5類函數(shù)：

具體函數(shù)的命名規(guī)則是：

函數(shù)名＝分布類型名稱+函數(shù)類型名稱(pdf、cdf、inv、stat、rnd)

其中，分布類型名稱如下：

分布類型MATLAB名稱

正態(tài)分布norm指數(shù)分布exp均勻分布unif

β分布beta

Γ分布gam對(duì)數(shù)正態(tài)分布lognrayleigh分布

raylweibull分布weib二項(xiàng)分布binoPoisson分布poiss幾何分布geo超幾何分布hyge離散均勻分布unid負(fù)二項(xiàng)分布nbin

例如，normpdf、normcdf、norminv、normstat和normrnd分別是正態(tài)分布的概率密度、累積分布、逆累積分布、數(shù)字特征和隨機(jī)數(shù)生成函數(shù)。

關(guān)于這5類函數(shù)的語(yǔ)法，請(qǐng)?jiān)斠娪嘘P(guān)書籍。

快捷的學(xué)習(xí)可借助MATLAB的系統(tǒng)幫助，通過指令doc獲得具體函數(shù)的詳細(xì)信息，語(yǔ)法是

doc<函數(shù)名>例6

到某服務(wù)機(jī)構(gòu)辦事總是要排隊(duì)等待的。設(shè)等待時(shí)間T是服從指數(shù)分布的隨機(jī)變量(單位：分鐘)，概率密度為

設(shè)某人一個(gè)月內(nèi)要到此辦事10次，若等待時(shí)間超過15分鐘，他就離去。求：

(1)恰好有兩次離去的概率；

(2)最多有兩次離去的概率；

(3)至少有兩次離去的概率；

(4)離去的次數(shù)占多數(shù)的概率。

解首先求任一次離去的概率，依題意

設(shè)10次中離去的次數(shù)為X，則。

>>p=1-expcdf(15,10)%任一次離去的概率p1=binopdf(2,10,p)%恰有兩次離去的概率q=binopdf([0:2],10,p);p2=sum(q)%最多有兩次離去的概率q=binopdf([0:1],10,p);p3=1-sum(q)%最少有兩次離去的概率q=binopdf([0:5],10,p);p4=1-sum(q)%離去的次數(shù)占多數(shù)的概率

p=0.2231p1=0.2972p2=0.6073p3=0.6899p4=0.0112例7某一急救中心在長(zhǎng)度為t的時(shí)間間隔內(nèi)收到的緊急呼救次數(shù)服從參數(shù)為t／2的泊松分布，而與時(shí)間間隔的起點(diǎn)無(wú)關(guān)(時(shí)間以小時(shí)計(jì))，求：

(1)在某一天中午12時(shí)至下午3時(shí)沒有收到緊急呼救的概率；

(2)某一天中午12時(shí)至下午5時(shí)至少收到1次緊急呼救的概率。

（1）>>P1=poisscdf(0,3/2)P1=0.2231或者>>P1=poisspdf(0,3/2)P1=0.2231中午12時(shí)到下午3時(shí)沒有收到緊急呼救的概率為0.2231。（2）>>P2=1-poisscdf(0,5/2)P2=0.9179中午12時(shí)至下午5時(shí)至少收到1次緊急呼救的概率為0.9179。解本題計(jì)算需調(diào)用函數(shù)poisscdf，其格式為poisscdf(x,λ)，返回的值。例8

某廠研發(fā)了一種新產(chǎn)品，現(xiàn)要設(shè)計(jì)它的包裝箱，要求每箱至少裝100件產(chǎn)品，且開箱驗(yàn)貨時(shí)，每箱至少裝有100件合格產(chǎn)品的概率不應(yīng)小于0.9，假設(shè)隨機(jī)裝箱時(shí)每箱中的不合格產(chǎn)品數(shù)服從參數(shù)為3的泊松分布。問：要設(shè)計(jì)的這種包裝箱，每箱至少應(yīng)裝多少件產(chǎn)品才能滿足要求？

解設(shè)每箱至少裝100+m件產(chǎn)品，X表示每箱中的不合格品數(shù)，則X服從參數(shù)為3的泊松分布，即，依題意，即要求按下面的不等式確定m。>>clear;clf,m=0;p=0;whilep<=0.9q=poisspdf([0:m],3);p=sum(q);m=m+1;endmm=6

計(jì)算結(jié)果表明當(dāng)m=6時(shí)，p>0.9。即設(shè)計(jì)的包裝箱每箱至少應(yīng)裝106件產(chǎn)品。

例9

某種重大疾病的醫(yī)療險(xiǎn)種，每份每年需交保險(xiǎn)費(fèi)100元，若在這一年中，投保人得了這種疾病，則每份可以得到索賠額10000元，假設(shè)該地區(qū)這種疾病的患病率為0.0002，現(xiàn)該險(xiǎn)種共有10000份保單，問：

(1)保險(xiǎn)公司虧本的概率是多少?

(2)保險(xiǎn)公司獲利不少于80萬(wàn)元的概率是多少?

解設(shè)表示這一年中發(fā)生索賠的份數(shù)，依題意，的統(tǒng)計(jì)規(guī)律可用二項(xiàng)分布來(lái)描述。由二項(xiàng)分布與泊松分布的近似計(jì)算關(guān)系有近似服從參數(shù)為2的泊松分布。當(dāng)索賠份數(shù)超過100份時(shí)，則保險(xiǎn)公司發(fā)生虧本，虧本的概率為當(dāng)索賠份數(shù)不超過20份時(shí)，則保險(xiǎn)公司獲利就不少于80萬(wàn)元，其概率為>>[p]=poisspdf([0:19],2);%計(jì)算出20個(gè)泊松分布概率值

或[p]=binopdf([0:19],10000,0.0002);%按二項(xiàng)分布計(jì)算

p2=sum(p)

%求出保險(xiǎn)公司獲利不少于80萬(wàn)元的概率

p2=1.0000>>[p]=poisspdf([0:100],2);%計(jì)算101個(gè)泊松分布概率值

或[p]=binopdf([0:100],10000,0.0002);%按二項(xiàng)分布計(jì)算

p1=1-sum(p)%求出保險(xiǎn)公司虧本的概率

p1=0.0000

例10

設(shè)，求

，。

本題計(jì)算正態(tài)分布的累積概率值，調(diào)用函數(shù)normcdf，其格式為normcdf(x,μ,σ)，返回的值。解：>>p1=normcdf(6,4,3)-normcdf(3,4,3)p1=0.3781>>p2=1-normcdf(3,4,3)p2=0.6306例11

繪制正態(tài)分布的密度函數(shù)、分布函數(shù)曲線，并求均值與方差。

解：>>clearmu=2.5;sigma=0.6;x=(mu-4*sigma):0.005:(mu+4*sigma);y=normpdf(x,mu,sigma);f=normcdf(x,mu,sigma);plot(x,y,'-g',x,f,':b')[M,V]=normstat(mu,sigma)legend('pdf','cdf',-1)M=2.5000V=0.3600

從圖中可以看出，正態(tài)密度曲線是關(guān)于x＝μ對(duì)稱的鐘形曲線(兩側(cè)在μ±σ處各有一個(gè)拐點(diǎn))，正態(tài)累積分布曲線當(dāng)x＝μ時(shí)F(x)＝0.5。例12

觀察正態(tài)分布參數(shù)對(duì)密度曲線的影響。解：>>clearmu1=2.5;mu2=3;sigma1=0.5;sigma2=0.6;x=(mu2-4*sigma2):0.01:(mu2+4*sigma2);y1=normpdf(x,mu1,sigma1);%考察均值的影響y2=normpdf(x,mu2,sigma1);y3=normpdf(x,mu1,sigma1);%考察方差的影響y4=normpdf(x,mu1,sigma2);subplot(1,2,1)%考察結(jié)果的可視化plot(x,y1,'-g',x,y2,'-b')xlabel('\fontsize{12}μ1<μ2,σ1=σ2')legend('μ1','μ2')subplot(1,2,2)plot(x,y3,'-g',x,y4,'-b')xlabel('\fontsize{12}μ1=μ2,σ1<σ2')legend('σ1','σ2')例13

正態(tài)分布參數(shù)μ和σ對(duì)變量x取值規(guī)律的約束——3σ準(zhǔn)則。解：>>clear,clf%（標(biāo)準(zhǔn)）正態(tài)分布密度曲線下的面積X=linspace(-5,5,100);Y=normpdf(X,0,1);yy=normpdf([-3,-2,-1,0,1,2,3],0,1);plot(X,Y,'k-',[0,0],[0,yy(4)],'c-.')holdonplot([-2,-2],[0,yy(2)],'m:',[2,2],[0,yy(6)],'m:',[-2,-0.5],[yy(6),yy(6)],'m:',[0.5,2],[yy(6),yy(6)],'m:')plot([-1,-1],[0,yy(3)],'g:',[1,1],[0,yy(5)],'g:',[-1,-0.5],[yy(5),yy(5)],'g:',[0.5,1],[yy(5),yy(5)],'g:')plot([-3,-3],[0,yy(1)],'b:',[3,3],[0,yy(7)],'b:',[-3,-0.5],[yy(7),yy(7)],'b:',[0.5,3],[yy(7),yy(7)],'b:')holdofftext(-0.5,yy(6)+0.005,'\fontsize{14}95.44%')text(-0.5,yy(5)+0.005,'\fontsize{14}68.26%')text(-0.5,yy(7)+0.005,'\fontsize{14}99.74%')text(-3.2,-0.03,'\fontsize{10}μ-3σ')text(-2.2,-0.03,'\fontsize{10}μ-2σ')text(-1.2,-0.03,'\fontsize{10}μ-σ')text(-0.05,-0.03,'\fontsize{10}μ')text(0.8,-0.03,'\fontsize{10}μ+σ')text(1.8,-0.03,'\fontsize{10}μ+2σ')text(2.8,-0.03,'\fontsize{10}μ+3σ')例14

標(biāo)準(zhǔn)正態(tài)分布α分位數(shù)的概念圖示。解>>%α分位數(shù)示意圖（標(biāo)準(zhǔn)正態(tài)分布，α=0.05）clear,clfdata=normrnd(0,1,300,1);xalpha1=norminv(0.05,0,1);xalpha2=norminv(0.95,0,1);xalpha3=norminv(0.025,0,1);xalpha4=norminv(0.975,0,1);subplot(3,1,1)capaplot(data,[-inf,xalpha1]);axis([-3,3,0,0.45])subplot(3,1,2)capaplot(data,[xalpha2,inf]);axis([-3,3,0,0.45])subplot(3,1,3)capaplot(data,[-inf,xalpha3]);axis([-3,3,0,0.45])holdoncapaplot(data,[xalpha4,inf]);axis([-3,3,0,0.45])holdoffxalpha1

xalpha2

xalpha3

xalpha4xalpha1=-1.6449xalpha2=1.6449xalpha3=-1.9600xalpha4=1.9600數(shù)理統(tǒng)計(jì)基礎(chǔ)Matlab統(tǒng)計(jì)工具箱中常見的統(tǒng)計(jì)命令1、基本統(tǒng)計(jì)量對(duì)于隨機(jī)變量x，計(jì)算其基本統(tǒng)計(jì)量的命令如下：均值：mean(x)標(biāo)準(zhǔn)差：std(x)中位數(shù)：median(x)方差：var(x)偏度：skewness(x)峰度：kurtosis(x)2、頻數(shù)直方圖的描繪A、給出數(shù)組data的頻數(shù)表的命令為：[N,X]=hist(data,k)

此命令將區(qū)間[min(data),max(data)]分為k個(gè)小區(qū)間(缺省為10)，返回?cái)?shù)組data落在每一個(gè)小區(qū)間的頻數(shù)N和每一個(gè)小區(qū)間的中點(diǎn)X。B、描繪數(shù)組data的頻數(shù)直方圖的命令為：hist(data,k)3、參數(shù)估計(jì)A、對(duì)于正態(tài)總體，點(diǎn)估計(jì)和區(qū)間估計(jì)可同時(shí)由以下命令獲得：[muhat,sigmahat,muci,sigmaci]=normfit(x,alpha)此命令在顯著性水平alpha下估計(jì)x的參數(shù)（alpha缺省值為5%），返回值muhat是均值的點(diǎn)估計(jì)值，sigmahat是標(biāo)準(zhǔn)差的點(diǎn)估計(jì)值，muci是均值的區(qū)間估計(jì)，sigmaci是標(biāo)準(zhǔn)差的區(qū)間估計(jì)。B、對(duì)其他分布總體，兩種處理辦法：一是取容量充分大的樣本，按中心極限定理，它近似服從正態(tài)分布，仍可用上面估計(jì)公式計(jì)算；二是使用特定分布總體的估計(jì)命令，常用的命令如：[muhat,muci]=expfit(x,alpha)[lambdahat,lambdaci]=poissfit(x,alpha)[phat,pci]=weibfit(x,alpha)4、正態(tài)總體假設(shè)檢驗(yàn)A、單總體均值的z檢驗(yàn)：

[h,sig,ci]=ztest(x,m,sigma,alpha,tail)檢驗(yàn)數(shù)據(jù)x關(guān)于總體均值的某一假設(shè)是否成立，其中sigma為已知方差，alpha為顯著性水平，究竟檢驗(yàn)什么假設(shè)取決于tail的取值：tail=0，檢驗(yàn)假設(shè)“x的均值等于m”tail=1，檢驗(yàn)假設(shè)“x的均值大于m”tail=-1，檢驗(yàn)假設(shè)“x的均值小于m”tail的缺省值為0，alpha的缺省值為5%。返回值h為一個(gè)布爾值，h=1表示可拒絕原假設(shè)，h=0表示不可拒絕原假設(shè)，sig為假設(shè)成立的概率，ci為均值的1-alpha置信區(qū)間。B、單總體均值的t檢驗(yàn)：

[h,sig,ci]=ttest(x,m,alpha,tail)C、雙總體均值的t檢驗(yàn)：

[h,sig,ci]=ttest2(x,y,alpha,tail)5、非參數(shù)檢驗(yàn)：總體分布的檢驗(yàn)Matlab統(tǒng)計(jì)工具箱提供了兩個(gè)對(duì)總體分布進(jìn)行檢驗(yàn)的命令：A、h=normplot(x)此命令顯示數(shù)據(jù)矩陣x的正態(tài)概率圖，如果數(shù)據(jù)來(lái)自于正態(tài)分布，則圖形顯示出直線形態(tài)，而其他概率分布函數(shù)顯示出曲線形態(tài)。B、h=weibplot(x)此命令顯示數(shù)據(jù)矩陣x的Weibull概率圖，如果數(shù)據(jù)來(lái)自于Weibull分布，則圖形顯示出直線形態(tài)，而其他概率分布函數(shù)顯示出曲線形態(tài)。例15

一道工序用自動(dòng)化車床連續(xù)加工某種零件，由于刀具損壞等會(huì)出現(xiàn)故障。故障是完全隨機(jī)的，并假定生產(chǎn)任一零件時(shí)出現(xiàn)故障機(jī)會(huì)均相同，工作人員是通過檢查零件來(lái)確定工序是否出現(xiàn)故障的?，F(xiàn)積累有100次故障紀(jì)錄，故障出現(xiàn)時(shí)該刀具完成的零件數(shù)如下：459362624542509584433748815505

612452434982640742565706593680926653164487734608428115359384452755251378147438882453886265977585975549697515628954771609402960885610292837473677358638699634555570844166061062484120447654564339280246687539790581621724531512577496468499544645764558378765666763217715310851

試觀察該刀具出現(xiàn)故障時(shí)完成的零件數(shù)屬于哪種分布？>>%數(shù)據(jù)輸入x1=[459362624542509584433748815505];x2=[612452434982640742565706593680];x3=[9266531644877346084281153593844];x4=[527552513781474388824538862659];x5=[77585975549697515628954771609];x6=[402960885610292837473677358638];x7=[699634555570844166061062484120];x8=[447654564339280246687539790581];x9=[621724531512577496468499544645];x10=[764558378765666763217715310851];x=[x1x2x3x4x5x6x7x8x9x10];%作頻數(shù)直方圖hist(x,10)[N,X]=hist(x,10)%分布的正態(tài)性檢驗(yàn)normplot(x)N=33714242214832X=1.0e+003*0.10420.21460.32500.43540.54580.65620.76660.87700.98741.0978>>%參數(shù)估計(jì)[muhat,sigmahat,muci,sigmaci]=normfit(x)muhat=594sigmahat=204.1301muci=553.4962

634.5038sigmaci=179.2276

237.1329刀具壽命服從正態(tài)分布，均值估計(jì)值為594，方差估計(jì)值為204.1301，均值的95%置信區(qū)間為[553.4962，634.5038]，方差的95%置信區(qū)間為[179.2276，237.1329]>>%假設(shè)檢驗(yàn)[h,sig,ci]=ttest(x,594)%已知刀具壽命服從正態(tài)分布，方差未知的情況下，檢驗(yàn)壽命均值是否等于594。h=0sig=1ci=553.4962

634.5038檢驗(yàn)結(jié)果：布爾變量h=0，表示不可拒絕原假設(shè)，說明假設(shè)壽命均值等于594是合理的。

95%置信區(qū)間為[553.4962，634.5038]完全包括594，估計(jì)精度較高。sig=1遠(yuǎn)超過0.05，不可拒絕原假設(shè)所以可以認(rèn)為刀具平均壽命為594（件）例16用模擬試驗(yàn)的方法直觀地驗(yàn)證教材§6.2抽樣分布定理一的結(jié)論。

假定變量，用隨機(jī)數(shù)生成的方法模擬對(duì)的500次簡(jiǎn)單隨機(jī)抽樣，每個(gè)樣本的容量為16。利用這500×16個(gè)樣本數(shù)據(jù)直觀地驗(yàn)證樣本均值的抽樣分布為均值等于60、方差等于25／16的正態(tài)分布，即解>>%1、用隨機(jī)數(shù)生成的方法模擬簡(jiǎn)單隨機(jī)抽樣x=[];%生成一個(gè)存放樣本數(shù)據(jù)的空表（維數(shù)可變的動(dòng)態(tài)矩陣）forbyk=1:500%循環(huán)控制，循環(huán)執(zhí)行下面的指令500次，本例中相當(dāng)于500次抽樣

xx=normrnd(60,5,16,1);%生成一個(gè)來(lái)自N(60,25)的容量為16的樣本（列向量）

x=[x,xx];%將樣本數(shù)據(jù)逐列存入數(shù)表x，可從matlab的變量瀏覽器（workspace）中觀察這個(gè)數(shù)表end

%2、計(jì)算每個(gè)樣本的樣本均值（1~500）xmean=mean(x);%可從變量瀏覽器中觀察這500個(gè)數(shù)據(jù)

%3、繪制500個(gè)樣本均值數(shù)據(jù)的直方圖k=ceil(1.87*(length(x)-1)^(2/5));%確定分組數(shù)h=histfit(xmean,k);%繪制附正態(tài)參考曲線的數(shù)據(jù)直方圖set(h(1),'FaceColor','c','EdgeColor','w')%修飾，設(shè)置直方圖線條顏色與填充色

%4、用這500個(gè)樣本均值數(shù)據(jù)驗(yàn)證樣本均值的均值和方差M=mean(xmean)%求（1~500）樣本的樣本均值的均值V=var(xmean)%求（1~500）樣本的樣本均值的方差M=59.9879V=1.4129M=60.0117V=1.3900M=59.9749V=1.5158M=59.9929V=1.5757M=59.8809V=1.6855…………例17觀察：用binornd模擬5000次投球過程，觀察小球堆積的情況。>>clear;clf,n=5;p=0.5;m=5000;x=[0:1:n]rand('seed',3)R=binornd(n,p,1,m);%模擬服從二項(xiàng)分布的隨機(jī)數(shù)，相當(dāng)于模擬投球m次forI=1:n+1%開始計(jì)數(shù)

k=[];k=find(R==(I-1));%find是一個(gè)有用的指令，本語(yǔ)句的作用是找出R中等于(I-1)元素下標(biāo)，并賦予向量k中

h(I)=length(k)/m;%計(jì)算落于編號(hào)(I-1)的格子中的小球頻率endbar(x,h),axis([-1601])%畫頻率圖title('\fontsize{18}\fontname{華文新魏}5000次投球小球堆積的頻率圖')>>f=binopdf(x,n,p),bar(x,f),axis([-1601])title('\fontsize{18}\fontname{華文新魏}B(5,0.5)理論分布圖')例18

利用隨機(jī)數(shù)樣本驗(yàn)證中心極限定理。

獨(dú)立同分布的隨機(jī)變量的和的極限分布服從正態(tài)分布，通過產(chǎn)生容量為n的poiss分布和exp分布的樣本，研究其和的漸近分布。

算法如下：①產(chǎn)生容量為n的獨(dú)立同分布的隨機(jī)數(shù)樣本，得其均值和標(biāo)準(zhǔn)差；②將隨機(jī)數(shù)樣本和標(biāo)準(zhǔn)化；③重復(fù)①、②；④驗(yàn)讓所得標(biāo)準(zhǔn)化的隨機(jī)數(shù)樣本和是否服從標(biāo)準(zhǔn)正態(tài)分布。

具體程序如下：>>clearn=2000;means=0;s=0;y=[];lamda=4;a=lamda;fori=1:nr=poissrnd(a,n,1);%可換成r=exprnd(a,n,1)；

means=mean(r);%計(jì)算樣本均值

s=std(r);%計(jì)算樣本標(biāo)準(zhǔn)差

y(i)=sqrt(n).*(means-a)./sqrt(s);endnormplot(y);%分布的正態(tài)性檢驗(yàn)title('poiss分布，中心極限定理')例19在同一坐標(biāo)軸上畫box圖，并對(duì)兩個(gè)班的成績(jī)進(jìn)行初步的分析比較。

兩個(gè)教學(xué)班各30名同學(xué)，在數(shù)學(xué)課程上，A班用新教學(xué)方法組織教學(xué)，B班用傳統(tǒng)方法組織教學(xué)，現(xiàn)得期末考試成績(jī)?nèi)缦隆?/p>

A：82，92，77，62，70，36，80，100，74，64，63，56，72，78，68，65，72．70，58，92，79,92，65，56，85，73，61，71，42，89

B：57，67，64，54，77，65，71，58，59，69，67，84，63，95，81，46，49,60,64，66，74，55，58，63，65，68，76，72，48，72解>>clear

x=[82,92,77,62,70,36,80,100,74,64,63,56,72,78,68,65,72,70,58,92,79,92,65,56,85,73,61,71,42,89;57,67,64,54,77,65,71,58,59,69,67,84,63,95,81,46,49,60,64,66,74,55,58,63,65,68,76,72,48,72];

boxplot(x')

從圖中直觀地看出，兩個(gè)班成績(jī)的分布是正態(tài)（對(duì)稱）的，A班成績(jī)較為分散(方差大)，B班成績(jī)則較集中(方差小)。A班成績(jī)明顯高于B班(均值比較．并且A班25％低分段上限接近B班中值線，A班中值線接近B班25％高分段下限)。A班的平均成績(jī)約為70分(中值)，B班約為65分(中值)。A班有一名同學(xué)的成績(jī)過低(離群)，而B班成績(jī)優(yōu)秀的只有一人(離群)。

需要注意的是，從圖中我們不能得出新教學(xué)方法一定優(yōu)于傳統(tǒng)教學(xué)方法的結(jié)論，因?yàn)槲覀儾⒉恢纼蓚€(gè)班級(jí)原有的數(shù)學(xué)基礎(chǔ)是怎樣的。三、MATLAB也為常用的三大統(tǒng)計(jì)分布提供了相應(yīng)的pdf、cdf、inv、stat、rnd類函數(shù)，具體分布類型函數(shù)名稱如下：

分布類型MATLAB名稱

分布chi2t分布tF分布f非中心分布ncx2非中心t分布

nct非中心F分布ncf例20分布的密度函數(shù)曲線。解：

>>%繪制不同自由度的卡方分布概率密度曲線clear,clfX=linspace(0,20,100);Y1=chi2pdf(X,1);%自由度等于1Y2=chi2pdf(X,3);%自由度等于3Y3=chi2pdf(X,6);%自由度等于6plot(X,Y1,'-g',X,Y2,'-b',X,Y3,'-k')title('\fontsize{18}\fontname{華文新魏}不同自由度的{\chi}^2分布概率密度曲線的比較')text(0.6,0.6,'\fontsize{12}df:n=1')text(2.6