楊曉霞

（北京林業(yè)大學 理學院，中國 北京100083）

概率統(tǒng)計是應用性很強的學科,在教學中不僅要重視理論知識的傳授，同時應引導學生關注計算問題，熟悉統(tǒng)計軟件。為了達到這一目的，利用統(tǒng)計軟件完成一些任務，能夠激發(fā)學生的成就感，進而提高學習興趣。因此，結合教學內(nèi)容進行一些概率統(tǒng)計小實驗，讓學生熟悉統(tǒng)計軟件及常用的計算方法，是改進概率統(tǒng)計教學的一個重要部分。

本文選取在概率統(tǒng)計教學中遇到的三個典型問題，利用蒙特卡羅(Monte Carlo)方法，借助于R 軟件，進行詳細分析和模擬計算，以期對概率統(tǒng)計的實驗教學提供參考。

蒙特卡羅方法，是一種基于隨機數(shù)的統(tǒng)計模擬方法。 在一些統(tǒng)計問題中，有時很難給出完美的理論結果，這時可以通過隨機模擬來給出近似的結果。 一般需要根據(jù)已知的概率分布構造出相應的模擬樣本，用它們的樣本頻率代替相應的概率作統(tǒng)計分析和推斷[1]。隨著概率統(tǒng)計學的發(fā)展，蒙特卡羅方法已成為普遍應用的方法之一。 甚至有時即使得到了問題的精確結果，也需要通過隨機模擬來進行驗證。

R 軟件是一個有著強大的統(tǒng)計分析及作圖功能的軟件系統(tǒng)，它在GNU 協(xié)議（General Public License）下免費發(fā)行，由R 核心小組成員開發(fā)及維護[2]。由于其免費共享且功能強大的特點，R 軟件已經(jīng)成為科研工作者最常使用的統(tǒng)計軟件之一。

1 民航送客問題

例1 一輛民航送客車載有20 位乘客從機場開出，沿途有10 個車站可以下車。 到達一個車站后若無乘客下車就不停車。 X 表示停車的次數(shù)，求EX。（假設每個乘客在各個車站下車是等可能的，且各乘客是否下車相互獨立）[3]。

而由分析可知P（Xi=0）=0.920，P（Xi=1）=1-0.920，i=1，2，…,10. 易計算得EXi=1-0.920， 所以EX=EX1+EX2+…+EX10=10×（1-0.920）=8.784（次）

分析：此題屬于數(shù)學期望性質的應用的題目，解法非常巧妙。但是在教學過程中發(fā)現(xiàn)， 學生們總試圖直接求出X 的分布， 然后計算期望。 然而X 的分布并不容易得到，因此解題過程進行不下去，使得學生的心里感覺不踏實。 而利用蒙特卡羅方法，我們很容易就可以模擬出X 的分布，進而得到EX 的估計值。

模擬算法步驟：

（1）從1～10 中有放回隨機抽取20 個數(shù)；（20 個乘客各自下車的車站數(shù)）

（2）若在20 個數(shù)中有m（0≤m≤10）個不相同的數(shù)，則X=m；（m個車站有人下車）

（3）將步驟(1)、(2) 重復n=10000 次，得到X 的10000 個值；

（4）求10000 個X 值的平均值，即為EX 的估計值。

R 軟件的模擬程序：

station=seq(1,10)

sp=rep(1/10,10)

n=10000

x=rep(0,n)

for ( i in 1:n){

x_sample=sample(station, 20, replace=TRUE, sp)

x[i]=sum(station %in% x_sample)

}

table(x)/n

mean(x)

表1 是停車次數(shù)X 的模擬分布。 這里， 雖然X 還可能取值1，2，3，4，但實際上相應的概率非常小，接近于0，因此沒有在表1 中顯示出來。

表1 重復n＝10000 次時停車次數(shù)X 的模擬分布

從表2 可看出模擬計算結果與理論值很接近，且重復次數(shù)n 越大與理論值越接近。

表2 重復n次時EX 的模擬結果

2 估計量的有效性問題

例2 已知隨機變量X 服從參數(shù)為λ 的泊松分布，由簡單隨機樣本X1，X2，…，Xn給出未知參數(shù)λ 的無偏估計量[4]。

模擬算法步驟：

（1）設定λ=λ0（=0.5），從參數(shù)為λ0的Poisson 分布中隨機抽取m（=100）個隨機數(shù)；

利用R 軟件中的函數(shù)rpois(),mean(),var(),sd()等同樣也可以寫出模擬程序來。

從結果看，兩個估計量的均值都與λ0=0.5 非常接近，但是S2的標準差要比的標準差大。 我們進一步對不同的參數(shù)λ 重復上述過程，結果如表3。 可看出，隨著參數(shù)λ 的增大，兩個估計量的標準差都在增大，但S2的標準差總比的標準差大，所以要比S2更有效。

?

3 用中心極限定理估算概率的問題

例3 一顆骰子獨立地擲四次， 出現(xiàn)的點數(shù)總和記為X， 估計P（10＜X≤18）=_________。

解：令Yi表示第i 次擲出骰子的點數(shù)，i=1，…,4，則易得E（Yi）=

分析：在解此題的過程中，學生常會產(chǎn)生疑問，四個隨機變量相加是否個數(shù)太少？ 估計的概率值是否接近真實值呢？ 這也正是蒙特卡羅方法可以解答的問題[5]。這里，每顆骰子擲得的點數(shù)服從離散型均勻分布：

模擬算法步驟：

（1）從上述離散型隨機變量分布中隨機抽取4 個隨機數(shù)，求和得到X；

（2）將步驟（1）重復n=10000 次，得到（X1，X2，…，X10000）；

（3）數(shù)出其中滿足10＜Xi≤18 的個數(shù)，除以n 即得到P（10＜X≤18）的估計值。

利用R 軟件中的函數(shù)sample(),sum(),length()等可以寫出模擬程序。 表4 給出了重復n 次時的模擬結果。

?

從表4 模擬的結果看，模擬概率值與由中心極限定理估計的概率0.796 相差不大，說明此例用中心極限定理進行估算是適用的。更進一步，雖然較繁瑣，但我們?nèi)菀椎贸鯴 的確切分布，并計算出P（10＜X≤18）的真實值約為0.7438。

在概率統(tǒng)計課程的學習中，通過學習蒙特卡羅方法，學生不僅能提出問題，也能通過隨機模擬自己解決問題。當然，這些只是起到一個引導的作用。 我們的目的是為學生提供一個有用的工具，并通過例子激發(fā)學生的興趣， 相信在需要的時候他們會自發(fā)學習和使用R 軟件中眾多的函數(shù)和功能。

［1］薛毅，陳立萍.統(tǒng)計建模與R 軟件[M].北京：清華大學出版社，2007，4：1-525.

［2］R Development Core Team(2005). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL:http://www.R-project.org[OL].

［3］盛驟,謝式千,潘承毅.概率論與數(shù)理統(tǒng)計[M].高等教育出版社,2008,6:1-362.

［4］賈乃光,張青,李永慈.數(shù)理統(tǒng)計[M].中國林業(yè)出版社,2006,7:1-215.

［5］何書元.概率論與數(shù)理統(tǒng)計[M].高等教育出版社,2006,6:1-372.