引言

而MATLAB軟件具有簡單易學、易操作和繪圖功能強等特點，利用MATLAB軟件的圖形可視功能將概率統計的內容用圖形表示出來，通過圖形讓學生加深理解，以達到事半功倍的效果。

概率論與數理統計知識比較抽象，邏輯性較強。因此，建議讓學生結合理論和公式推導，進行數值試驗和相關調查，直觀地感受數學概念和理論，從而提高學生解決實際問題的信心和能力。概率論

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

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

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

間，生成方式為正態分布。

3.randperm(m)：生成一個1~m的隨機整數排列。

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

5.取整函數系列：

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

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

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

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

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

7.prod(x)：向量x的所有分量元素的積。一、MATLAB常用的與隨機數產生相關的函數：示例：

>>rand(1)%生成一個(0,1)間的隨機數

ans=0.8147

>>rand(2,2)%生成一個2×2階(0,1)間的隨機數矩陣

ans=0.91340.0975

0.63240.2785

>>randperm(5)%生成一個1~5的隨機整數排列

ans=41523

>>a=[1242332];

unique(a)

ans=1234

例1

隨機投擲均勻硬幣，觀察國徽朝上與國徽

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

試驗次數n300050001萬2萬3萬國徽朝上頻率0.50400.50060.48790.49990.5046國徽朝下頻率0.49600.49940.51210.50010.4954試驗次數n5萬10萬100萬100萬1億國徽朝上頻率0.50210.49990.49990.50010.5000國徽朝下頻率0.49790.50010.50010.49990.5000可見當時，解記事件為第i個人拿到自已槍，事件為第i個人沒拿到自己槍，易知：又記為沒有一個人拿到自己槍的概率。

有乘法公式可知：例2

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

所以

特別地，當n較大時，。因此，可隨機模擬出沒有人拿到自己槍的頻率，根據頻率的穩定性，近似當做概率，然后去估計自然常數e。算法如下：

1、產生n個隨機數的隨機序列；

2、檢驗隨機列與自然列是否至少有一個配對；

3、對沒有一個配對的序列進行累積p；

4、重復1、2、3步m

次；

5、估計。

具體程序及相關結果為（注：自然常數

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模擬次數m400004000040000人數n100020005000e2.71552.70822.7202模擬次數m400040000400000人數n505050e2.73792.73132.7194

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

故可得的近似計算公式，其中n為隨機試驗次數，m為針與平行線相交的次數。例3Buffon投針實驗

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

>>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)試驗次數n5千1萬10萬100萬1000萬針長l/平行間距d3/103/103/103/103/10相交頻率0.18360.19710.18870.19050.1912π的近似值3.26803.04413.17983.14983.1387試驗次數n5千1萬10萬100萬1000萬針長l/平行間距d2/52/52/52/52/5相交頻率0.24960.25620.25490.25440.2543π的近似值3.20513.12263.13863.14513.1433試驗次數n5千1萬10萬100萬1000萬針長l/平行間距d1/21/21/21/21/2相交頻率0.32540.31480.31580.31780.3183π的近似值3.07313.17663.16673.14703.1417試驗次數n5千1萬10萬100萬1000萬針長l/平行間距d4/54/54/54/54/5相交頻率0.51420.51340.50860.50930.5093π的近似值3.11163.11653.14603.14183.1418試驗次數n5千1萬10萬100萬1000萬針長l/平行間距d17/2017/2017/2017/2017/20相交頻率0.54320.54520.54200.54120.5410π的近似值3.12963.11813.13663.14133.1426試驗次數n5千1萬10萬100萬1000萬針長l/平行間距d9/109/109/109/109/10相交頻率0.58600.57000.57560.57330.5731π的近似值3.07173.15793.12723.13953.1410例4在100個人的團體中，不考慮年齡差異，研究是否有兩個以上的人生日相同。假設每人的生日在一年365天中的任意一天是等可能的，那么隨機找n個人(不超過365人)。

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

(2)近似計算在30名學生的一個班中至少有兩個人生日相同的概率是多少？解:(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('人數'),ylabel('概率')legend('生日各不相同的概率','至少兩人生日相同的概率')axis([0100-0.11.199]),gridonp1(30)=0.7063,p1(60)=0.9941

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

隨機產生30個正整數，代表一個班30名學生的生日，然后觀察是否有兩人以上生日相同。當30個人中有兩人生日相同時，輸出“1”，否則輸出“0”。如此重復觀察100次，計算出這一事件發生的頻率。

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

y=0;x=1+fix(365*rand(1,30));%產生30個隨機數

fori=1:29%用二重循環尋找30個隨機數中是否有相同數

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

Galton釘板試驗是由英國生物統計學家和人類學家Galton設計的。故而得名。

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

12345678高爾頓釘板試驗小球最后落入的格數?記小球向右落下的次數為

則記小球向左落下的次數為

則符號函數,大于0返回1,小于0返回-1,等于0返回0

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

-6

-5

-4

-3

-2

-1

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

模擬Galton釘板試驗的步驟：

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

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

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

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

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

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

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

movien(n)：創建動畫矩陣；制作動畫矩陣數據；

Getframe：拷貝動畫矩陣；

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

M文件如下：解：>>clear,clf,m=100;n=5;y0=2;%設置參數ballnum=zeros(1,n+1);p=0.5;q=1-p;fori=n+1:-1:1%創建釘子的坐標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);%動畫開始，模擬小球下落路徑fori=1:ms=rand(1,n);%產生n個隨機數

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

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);%小球下落點的坐標

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

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

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

mm(i)=getframe;%存儲動畫數據

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

具體函數的命名規則是：

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

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

分布類型MATLAB名稱

正態分布norm指數分布exp均勻分布unif

β分布beta

Γ分布gam對數正態分布lognrayleigh分布

raylweibull分布weib二項分布binoPoisson分布poiss幾何分布geo超幾何分布hyge離散均勻分布unid負二項分布nbin

例如，normpdf、normcdf、norminv、normstat和normrnd分別是正態分布的概率密度、累積分布、逆累積分布、數字特征和隨機數生成函數。

關于這5類函數的語法，請詳見有關書籍。

快捷的學習可借助MATLAB的系統幫助，通過指令doc獲得具體函數的詳細信息，語法是

doc<函數名>例6

到某服務機構辦事總是要排隊等待的。設等待時間T是服從指數分布的隨機變量(單位：分鐘)，概率密度為

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

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

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

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

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

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

設10次中離去的次數為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)%離去的次數占多數的概率

p=0.2231p1=0.2972p2=0.6073p3=0.6899p4=0.0112例7某一急救中心在長度為t的時間間隔內收到的緊急呼救次數服從參數為t／2的泊松分布，而與時間間隔的起點無關(時間以小時計)，求：

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

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

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

某廠研發了一種新產品，現要設計它的包裝箱，要求每箱至少裝100件產品，且開箱驗貨時，每箱至少裝有100件合格產品的概率不應小于0.9，假設隨機裝箱時每箱中的不合格產品數服從參數為3的泊松分布。問：要設計的這種包裝箱，每箱至少應裝多少件產品才能滿足要求？

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

計算結果表明當m=6時，p>0.9。即設計的包裝箱每箱至少應裝106件產品。

例9

某種重大疾病的醫療險種，每份每年需交保險費100元，若在這一年中，投保人得了這種疾病，則每份可以得到索賠額10000元，假設該地區這種疾病的患病率為0.0002，現該險種共有10000份保單，問：

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

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

解設表示這一年中發生索賠的份數，依題意，的統計規律可用二項分布來描述。由二項分布與泊松分布的近似計算關系有近似服從參數為2的泊松分布。當索賠份數超過100份時，則保險公司發生虧本，虧本的概率為當索賠份數不超過20份時，則保險公司獲利就不少于80萬元，其概率為>>[p]=poisspdf([0:19],2);%計算出20個泊松分布概率值

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

p2=sum(p)

%求出保險公司獲利不少于80萬元的概率

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

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

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

p1=0.0000

例10

設，求

，。

本題計算正態分布的累積概率值，調用函數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

繪制正態分布的密度函數、分布函數曲線，并求均值與方差。

解：>>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

從圖中可以看出，正態密度曲線是關于x＝μ對稱的鐘形曲線(兩側在μ±σ處各有一個拐點)，正態累積分布曲線當x＝μ時F(x)＝0.5。例12

觀察正態分布參數對密度曲線的影響。解：>>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)%考察結果的可視化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

正態分布參數μ和σ對變量x取值規律的約束——3σ準則。解：>>clear,clf%（標準）正態分布密度曲線下的面積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

標準正態分布α分位數的概念圖示。解>>%α分位數示意圖（標準正態分布，α=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數理統計基礎Matlab統計工具箱中常見的統計命令1、基本統計量對于隨機變量x，計算其基本統計量的命令如下：均值：mean(x)標準差：std(x)中位數：median(x)方差：var(x)偏度：skewness(x)峰度：kurtosis(x)2、頻數直方圖的描繪A、給出數組data的頻數表的命令為：[N,X]=hist(data,k)

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

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

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

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

一道工序用自動化車床連續加工某種零件，由于刀具損壞等會出現故障。故障是完全隨機的，并假定生產任一零件時出現故障機會均相同，工作人員是通過檢查零件來確定工序是否出現故障的?，F積累有100次故障紀錄，故障出現時該刀具完成的零件數如下：459362624542509584433748815505

612452434982640742565706593680926653164487734608428115359384452755251378147438882453886265977585975549697515628954771609402960885610292837473677358638699634555570844166061062484120447654564339280246687539790581621724531512577496468499544645764558378765666763217715310851

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

634.5038sigmaci=179.2276

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

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

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

假定變量，用隨機數生成的方法模擬對的500次簡單隨機抽樣，每個樣本的容量為16。利用這500×16個樣本數據直觀地驗證樣本均值的抽樣分布為均值等于60、方差等于25／16的正態分布，即解>>%1、用隨機數生成的方法模擬簡單隨機抽樣x=[];%生成一個存放樣本數據的空表（維數可變的動態矩陣）forbyk=1:500%循環控制，循環執行下面的指令500次，本例中相當于500次抽樣

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

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

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

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

%4、用這500個樣本均值數據驗證樣本均值的均值和方差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);%模擬服從二項分布的隨機數，相當于模擬投球m次forI=1:n+1%開始計數

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

h(I)=length(k)/m;%計算落于編號(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

利用隨機數樣本驗證中心極限定理。

獨立同分布的隨機變量的和的極限分布服從正態分布，通過產生容量為n的poiss分布和exp分布的樣本，研究其和的漸近分布。

算法如下：①產生容量為n的獨立同分布的隨機數樣本，得其均值和標準差；②將隨機數樣本和標準化；③重復①、②；④驗讓所得標準化的隨機數樣本和是否服從標準正態分布。

具體程序如下：>>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);%計算樣本均值

s=std(r);%計算樣本標準差

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

兩個教學班各30名同學，在數學課程上，A班用新教學方法組織教學，B班用傳統方法組織教學，現得期末考試成績如下。

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')

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

需要注意的是，從圖中我們不能得出新教學方法一定優于傳統教學方法的結論，因為我們并不知道兩個班級原有的數學基礎是怎樣的。三、MATLAB也為常用的三大統計分布提供了相應的pdf、cdf、inv、stat、rnd類函數，具體分布類型函數名稱如下：

分布類型MATLAB名稱

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

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

>>%繪制不同自由度的卡方分布概率密度曲線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