面向計算思維的蒙特卡羅C語言程序設計案例探究

顧麗紅，丁淑妍

（中國石油大學（華東） 計算機與通信工程學院，山東 青島 266580）

0 引 言

計算思維[1]概念由Jeannette M. Wing教授于2006年首次完整提出，隨后引起了國內計算機教育界的重視。教育部計算機專業(yè)教學指導委員會明確提出將計算思維列入計算機專業(yè)人才的專業(yè)基本能力[2]。如何培養(yǎng)非計算機專業(yè)大學生的計算思維能力，這是擺在大學計算機基礎教育工作者面前的難題。目前國內高校的C語言程序設計課程普遍存在注重知識基礎而忽視能力基礎的問題，能力的培養(yǎng)則需要引入更多能引起學生興趣的程序設計案例實踐[3-5]。近幾年相關的教研論文越來越多，大多以數(shù)值計算類為主，面向跨計算學科的廣義計算思維[6]的案例設計實踐類論文還不多見。計算機模擬通常綜合了計算機科學家和領域科學家的思維方式，在課程教學實踐中引入這類程序設計案例，有利于培養(yǎng)學生的計算思維能力。

1 計算思維與計算機模擬

進入信息和數(shù)據(jù)時代，計算思維是每個人的基本技能，就像邏輯思維、數(shù)學思維、設計思維一樣，是人們認識問題和解決問題的重要思維方式之一。計算思維能力是人們有意識地使用計算機科學家的思想、方法、技術、工具、資源、環(huán)境進行思考和活動的能力。如何全方位培養(yǎng)這種能力是教育工作者面臨的挑戰(zhàn)。

計算機模擬通常是指用抽象模型來模擬特定系統(tǒng)的計算機程序。人們在研究、設計、構造復雜系統(tǒng)時，往往需要設計制造一個模型來進行各種試驗。傳統(tǒng)方法是建立數(shù)學模型或實物模型，前者缺乏直觀性也不便于試驗，后者比較直觀但不夠經(jīng)濟和方便。隨著計算機的出現(xiàn)和處理能力的不斷增強，利用計算機來模擬各種復雜系統(tǒng)，成為重要的科學研究手段。計算機模擬的一般過程如圖1所示。

其中，①數(shù)學建模與形式化。明確模擬工作的目標，確定模擬方案的評價準則。②模擬建模。根據(jù)問題特點和模擬要求從穩(wěn)定性、精度和性能3個維度設計合適的算法。③程序設計。把模擬模型用計算機能執(zhí)行的程序來描述。④模擬運行。分析運行結果是否合適，通過前幾步查找問題，不斷修正，直到結果滿意。⑤模擬實驗和結果處理。每一個方案，采用不同隨機數(shù)序列重復多次，采用數(shù)理統(tǒng)計方法分析模擬結果的統(tǒng)計特征；不同方案，采用相同隨機數(shù)序列進行模擬運行，消除因隨機數(shù)序列不同而引起的差異；模擬結果圖形化、可視化，乃至虛擬現(xiàn)實環(huán)境。

圖1 計算機模擬的一般處理過程

計算機模擬的設計原則包括：①分級模擬原則。采用與任務匹配的模型或算法，合理簡化問題，突出問題的關鍵。②效率兼顧原則。在模擬的不同階段，合理調整對準確度和速度的要求，提高模擬效率。③可信驗證原則。所有模擬結果應該有可信的驗證方法或依據(jù)。

在教學實踐中，讓學生充分理解計算機模擬過程中的算法構建、迭代修正、數(shù)理統(tǒng)計方法，以及結果圖形化和可視化的重要性，在提升學生興趣的實踐中培養(yǎng)能力。

2 蒙特卡羅方法思想

蒙特卡羅方法是一種隨機模擬方法，以概率和統(tǒng)計理論方法為基礎，使用隨機數(shù)（偽隨機數(shù)）解決計算問題，又稱為統(tǒng)計模擬法。因其具有概率統(tǒng)計特征而獲得數(shù)學家馮·諾依曼用賭城蒙特卡羅來命名。蒙特卡羅方法起源于1777年法國數(shù)學家蒲豐（Buffon）提出用投針試驗計算圓周率π，進而推廣到用概率法計算不規(guī)則圖形面積，甚至用小規(guī)模抽樣調查進行民意測驗，來預測競選的優(yōu)勝者，也是同樣的思想。當然科學計算中的問題比這復雜很多，比如金融領域的交易風險評估，問題的維數(shù)（隨機變量的個數(shù)）高達數(shù)百、幾千，問題的難度成指數(shù)增長，傳統(tǒng)數(shù)值方法難以勝任，而蒙特卡羅方法因計算復雜性不依賴于維數(shù)則可以很好地對付維數(shù)災難。

蒙特卡羅方法具有簡單、快速、節(jié)省存儲單元的優(yōu)點，省卻了一般人難以掌握的復雜的數(shù)學推導和演算過程，而且不受問題的幾何形狀復雜性影響，適用于處理大型復雜問題。隨著計算機處理復雜系統(tǒng)的能力日益增強，蒙特卡羅方法的應用也越來越廣泛。它不僅較好地解決了多重積分計算、微分方程求解、積分方程求解、特征值計算和非線性方程組求解等高難度和復雜的數(shù)學計算問題，而且在計算物理學、核物理、金融工程學、宏觀經(jīng)濟學、生物醫(yī)學、可靠性及計算機科學等廣泛的領域都得到成功的應用。

在教學實踐中，通過引入蒲豐投針問題和不規(guī)則面積計算問題，學生對蒙特卡羅方法思想有了濃厚的興趣，為之后的C語言算法設計思想的養(yǎng)成打下基礎。

3 蒙特卡羅模擬案例設計

蒙特卡羅模擬的基本思想：當所求問題的解是某個事件的概率，或是某個隨機變量的數(shù)學期望，或是與概率、數(shù)學期望有關的量時，通過某種試驗的方法，得出該事件發(fā)生的頻率或者該隨機變量的若干個具體觀察值的算術平均值，進而得到問題的解。采用蒙特卡羅方法進行計算機模擬的步驟：首先設計反映系統(tǒng)各部分運行時邏輯關系的邏輯框圖（模擬模型），然后通過具有各種概率分布的模擬隨機數(shù)來模擬隨機現(xiàn)象。

3.1 3種蒙特卡羅π的計算模擬

法國數(shù)學家蒲豐于18世紀提出投針問題：設有一個以平行且等距木紋鋪成的地板，現(xiàn)在隨意拋一支長度比木紋之間距離小的針，求針和其中一條木紋相交的概率，并以此概率可以近似計算圓周率。

蒲豐投針問題的數(shù)學建模：平行線距離為2a，針的長度為2l，且l≤a。設針投到地面上是隨機的，所以位置可以用二維隨機變量(x,θ)來描述，x為針中心的坐標，θ為針與平行線的夾角。任意投針，就是意味著x與θ都是任意取的，但x的范圍限于[0,a]，夾角θ的范圍限于[0, π]。在此情況下，針與平行線相交的數(shù)學條件是：x≤l×sinθ。x和θ均服從均勻分布，且相互獨立，通過(x,θ)的概率密度函數(shù)求出概率分布：P= 2l/πa，所以π的近似值為：2l/aP。

蒲豐投針C語言算法設計思路：

（1）隨機生成一個介于[0,a]表示針中心坐標的數(shù)x和一個介于[0, π]表示針與平行線的角度y；

（2）如果x≤l×siny，則表示針與平行線相交；

（3）計算π的近似值。

C語言程序設計實現(xiàn)如下：

蒲豐投針實驗的重要性在于它是第一個用幾何形式表達概率問題的例子，開創(chuàng)了使用隨機數(shù)處理確定性數(shù)學問題的先河。在教學實踐中，蒲豐實驗問題新穎奇妙，引起學生極大的興趣。為了達到舉一反三的效果，案例設計中引出蒙特卡羅思想計算π的另外兩種解法：圓中投點（幾何法）和隨機數(shù)互質法。

圓中投點問題描述：有一正方形木板，以及內切正方形內接圓盤，隨機向正方形中投點，求點落在圓盤中的概率，并以此概率近似計算圓周率。

圓中投點數(shù)學建模：有一個以（0, 0）為中心的邊長為2的正方形，以及這個正方形中半徑為1的內接圓，隨機向正方形中投點，求點落在內接圓中的概率。設點投到木板上是隨機的，所以位置可以用二維隨機變量(x,y)來描述。任意投點，就是意味著x與y都是任意取的，二者的范圍限于[0, 1]。在此情況下，落在內接圓中的數(shù)學條件是：x2+y2≤ 1。x和y均服從均勻分布，且相互獨立，求出概率分布：P= π/4，所以π的近似值為：4×P。

圓中投點C語言算法設計思路：

（1）隨機生成兩個介于[0, 1]區(qū)間的表示點坐標的x和y；

（2）如果x2+y2≤1，則表示點落在內接圓中；

（3）計算π的近似值。

C語言程序設計實現(xiàn)如下：

隨機數(shù)互質問題描述：R·查特在1904年發(fā)現(xiàn)，兩個隨意寫出的數(shù)中，互質的概率為6/π2，并以此概率近似計算圓周率。

隨機數(shù)互質數(shù)學建模：任意產(chǎn)生兩個整數(shù)x和y，x和y均服從均勻分布，且相互獨立，通過數(shù)學推導（推導過程較復雜）可以求出x和y互質的概率為：P= 6/π2，所以π的近似值為6/P的算術平方根。

隨機數(shù)互質C語言算法設計思路：

（1）隨機生成兩個介于[0, 32 767)的x和y；

（2）如果x和y互為質數(shù)，則統(tǒng)計變量累加1（注意特別處理生成數(shù)0的情況）；

（3）計算π的近似值。

C語言程序設計實現(xiàn)如下：

測試實驗的硬件環(huán)境：Intel i5 M540處理器，主頻2.53GHz，雙核四線程，3MB高速緩存。軟件環(huán)境：Window7旗艦版，CodeBlock 13.12 IDE環(huán)境，TDM-GCC v4.8.1編譯器。為了盡可能減少其他應用程序對測試的干擾，測試前盡可能關閉開機啟動項的應用軟件，比如殺毒軟件、下載工具等，并禁用網(wǎng)絡。采用多次測試取均值的方法可以降低測試誤差，本實驗取10次測試的平均值，實驗結果如表1。

眾所周知，圓周率π的真值在3.141 592 6和3.141 592 7之間，假設以中間值3.141 592 65為真值，我們發(fā)現(xiàn)當模擬的次數(shù)N為千萬次和億次的誤差分別為： 0.438‰和0.048‰， 0.235‰和0.023‰， 0.071‰和0.030‰，π的精度有了明顯的提高；當N為10億次時，精度沒有明顯提升，求更高精度的π需要尋求其他方法。蒙特卡羅方法的收斂速度與N的算術平方根成比例，意味著要使結果精度提高一位，應該增加一百倍的模擬計算工作量。

3.2 4種概率分布的隨機變量模擬

蒙特卡羅方法需要模擬隨機事件，即需要隨機數(shù)，那么計算機是如何產(chǎn)生隨機數(shù)的呢？隨機數(shù)是來自統(tǒng)計學的概念，真正的隨機數(shù)（真隨機數(shù)）是使用物理現(xiàn)象產(chǎn)生的，比如常見的彩票的搖號、擲錢幣、骰子、轉輪等，密碼學領域需要真隨機數(shù)。目前計算機產(chǎn)生的隨機數(shù)，是通過一個固定的、可以重復的算法產(chǎn)生，不是真正的隨機，但是具有類似于隨機數(shù)的統(tǒng)計特征，通常稱為偽隨機數(shù)，可以滿足蒙特卡羅方法的要求，所以默認計算機產(chǎn)生的隨機數(shù)是隨機的。

用蒙特卡羅方法模擬系統(tǒng)或過程時，需要不同概率分布的隨機數(shù)，通常先生成均勻分布的（0, 1）隨機數(shù)，這是其他分布隨機數(shù)的基礎。理論上說，具有連續(xù)分布的隨機數(shù)，通過數(shù)學變換或近似等方法，可以生成其他任意分布的隨機數(shù)。設μi(i=1, 2, …)是區(qū)間[0, 1]內的均勻分布的獨立隨機變量，而另一給定分布函數(shù)F(x)的隨機變量為則這一隨機變量xi可以由其反分布函數(shù)求得其抽樣值，即

均勻分布隨機變量算法設計思路：已知連續(xù)型隨機變量在有限區(qū)間內取值，則其概率密度函數(shù)為其分布函數(shù)為知所以

指數(shù)分布隨機變量算法設計思路：已知連續(xù)型隨機變量在有限區(qū)間內取值，則其概率密度函數(shù)為故其分布函數(shù)為知若(1-μ)是（0, 1）均勻分布隨機數(shù)，則可用下式簡化為

正態(tài)分布隨機變量算法設計思路：已知連續(xù)型隨機變量服從正態(tài)分布N(μ, σ2)，其概率密度函數(shù)為取區(qū)間（0, 1）上兩個均勻分布隨機數(shù)μ1和μ2，利用二元函數(shù)變換可得到個獨立的標準正態(tài)分布隨機變量。

表1 三種蒙特卡羅π的計算模擬結果

泊松分布隨機變量算法設計思路：泊松分布源自二項分布，在二項分布的伯努利試驗中，如果試驗次數(shù)n很大，二項分布的概率p很小，且乘積λ=np比較適中，則事件出現(xiàn)的次數(shù)的概率可以用泊松分布來逼近。泊松分布是離散型概率分布，表示固定尺度的連續(xù)區(qū)間上給定的事件發(fā)生次數(shù)的概率，其中n可以看作無窮大，泊松分布單位時間內隨機事件發(fā)生的次數(shù)滿足泊松分布，比如電話交換機接收到呼叫的次數(shù)、汽車站臺的候客人數(shù)、機器出現(xiàn)故障的次數(shù)、自然災害發(fā)生的次數(shù)、DNA序列的變異數(shù)等。

在C語言程序設計中，利用計算機產(chǎn)生的均勻分布隨機數(shù)，通過概率積分變換算法可以比較容易得到滿足其他概率分布的隨機變量。C語言中產(chǎn)生隨機數(shù)，通常用到兩個函數(shù)：srand()和rand()。srand()用來為計算機產(chǎn)生隨機數(shù)設置seed種子，否則每次程序運行產(chǎn)生的隨機數(shù)序列是一樣的。避免這種情況的通常方法是采用計算機的當前時間作為seed種子，因為時間是在不斷變化的。rand()函數(shù)隨機生成滿足均勻分布的[0,RAND_MAX]（RAND_MAX 為 0x7fff，即 32 767）之間的整數(shù)。如果要使范圍大一點，可以通過產(chǎn)生幾個隨機數(shù)的線性組合來實現(xiàn)任意范圍內的均勻分布隨機數(shù)，而且經(jīng)過一定的四則運算和取模運算，可以比較容易地得到任意區(qū)間的隨機變量，比如：均勻分布、指數(shù)分布、正態(tài)分布和泊松分布。

產(chǎn)生高精度均勻分布隨機數(shù)。用rand()函數(shù)產(chǎn)生隨機數(shù)，需要對精度有所認知。rand()的隨機數(shù)分辨率為32 767，兩個也就是65 534，如果需要要產(chǎn)生（-1 000, 1 000）之間且精度為4位小數(shù)點的均勻分布，則分辨率要求為1 000×10 000×2=20 000 000，這樣顯然遠遠不夠，但可以用兩個隨機數(shù)的乘法來達到精度要求，C語言程序設計片段如下。

利用[0, 1]區(qū)間的均勻分布可以產(chǎn)生指數(shù)分布隨機數(shù)，C語言程序設計片段如下。

利用[0, 1]區(qū)間的均勻分布可以產(chǎn)生正態(tài)分布隨機數(shù)（Box Muller方法），C語言程序設計片段如下。

利用[0, 1]區(qū)間的均勻分布可以產(chǎn)生泊松分布隨機數(shù)，Knuth算法設計思路：

1）l賦初值e-λ，k賦初值0，p賦初值1.0；

2）生成(0, 1)區(qū)間均勻分布隨機數(shù)u，p賦值為p*u，k累加1；

3）如果p>l，則重復2），直至不滿足條件；

4）返回k-1。

C語言程序設計片段如下：

實驗軟硬件環(huán)境同上，實驗結果如圖2（頻率曲線圖）。

從4種分布的頻率曲線圖可以看出，通過線性組合和四則運算，計算機產(chǎn)生的偽隨機數(shù)可以滿足蒙特卡羅方法的要求，在實際生產(chǎn)、生活和工程項目中可以大量應用。

4 結 語

蒙特卡羅模擬有很多專門的工具軟件，比如Mathmetica、MatLab、SPSS、CrystalBall等，其功能強大且容易使用。本文僅從C程序設計課程教學和實踐環(huán)節(jié)的需要，提出通過計算機模擬隨機事件引導學生對計算思維的認知和培養(yǎng)學生動手實踐的能力。大數(shù)據(jù)時代需要培養(yǎng)更多理解概率論和數(shù)理統(tǒng)計原理的融合型復合人才，建議大學計算機基礎課程教學工作者根據(jù)學生的專業(yè)背景設計更多跨學科領域且體現(xiàn)計算思維思想的案例庫進行交流和分享。下一步將在本文工作的基礎上，從大數(shù)據(jù)分析的角度，研究設計諸如投資風險分析、隨機庫存預測等進一步培養(yǎng)計算思維理能力的綜合型教學案例。

圖2 計算機生成的4種分布隨機數(shù)的頻率曲線圖

[1]Wing J M. Computational thinking[J]. Communications of the ACM, 2006, 49(3): 33-35.

[2]教育部高等學校計算機科學與技術教學指導委員會. 高等學校計算機科學與技術專業(yè)人才專業(yè)能力構成與培養(yǎng)[M]. 北京: 機械工業(yè)出版社, 2010.

[3]顧麗紅, 李傳秀, 吳少剛. 培養(yǎng)計算思維能力的矩陣乘法C語言程序設計案例探究[J]. 計算機教育, 2016(1): 149-152.

[4]姚天昉. 在程序設計課程中引入“計算思維”的實踐[J]. 中國大學教學, 2012(2): 61-62.

[5]陳文智, 陳越, 莊越挺. 面向系統(tǒng)設計能力培養(yǎng)的教學改革探索[J]. 計算機教育, 2013(20): 70-76.

[6]蔣宗禮. 計算思維之我見[J]. 中國大學教學, 2013(9): 5-10.