解析Monte-Carlo演算法
限時乾貨下載：
來源：http://www.cnblogs.com/leoo2sk/archive/2009/05/29/1491526.html
作者：張洋
1
引言
最近在和同學討論研究Six Sigma（六西格瑪）軟體開發方法及CMMI相關問題時，遇到了需要使用Monte-Carlo演算法模擬分布未知的多元一次概率密度分布問題。於是花了幾天時間，通過查詢相關文獻資料，深入研究了一下Monte-Carlo演算法，並以實際應用為背景進行了一些實驗。
在研究和實驗過程中，發現Monte-Carlo演算法是一個非常有用的演算法，在許多實際問題中，都有用武之地。目前，這個演算法已經在金融學、經濟學、工程學、物理學、計算科學及計算機科學等多個領域廣泛應用。而且這個演算法本身並不複雜，只要掌握概率論及數理統計的基本知識，就可以學會並加以應用。由於這種演算法與傳統的確定性演算法在解決問題的思路方面截然不同，作為計算機科學與技術相關人員以及程序員，掌握此演算法，可以開闊思維，為解決問題增加一條新的思路。
基於以上原因，我有了寫這篇文章的打算，一是回顧總結這幾天的研究和實驗，加深印象，二是和朋友們分享此演算法以及我的一些經驗。
這篇文章將首先從直觀的角度，介紹Monte-Carlo演算法，然後介紹演算法基本原理及數理基礎，最後將會和大家分享幾個基於Monte-Carlo方法的有意思的實驗。所以程序將使用C#實現。
閱讀本文需要有一些概率論、數理統計、微積分和計算複雜性的基本知識，不過不用太擔心，我將盡量避免過多的數學描述，並在適當的地方對於用到的數學知識進行簡要的說明。
2
Monte-Carlo演算法引導
首先，我們來看一個有意思的問題：在一個1平方米的正方形木板上，隨意畫一個圈，求這個圈的面積。
我們知道，如果圓圈是標準的，我們可以通過測量半徑r，然後用 S = pi * r^2 來求出面積。可是，我們畫的圈一般是不標準的，有時還特別不規則，如下圖是我畫的巨難看的圓圈。
圖1、不規則圓圈
顯然，這個圖形不太可能有面積公式可以套用，也不太可能用解析的方法給出準確解。不過，我們可以用如下方法求這個圖形的面積：
假設我手裡有一支飛鏢，我將飛鏢擲向木板。並且，我們假定每一次都能擲在木板上，不會偏出木板，但每一次擲在木板的什麼地方，是完全隨機的。即，每一次擲飛鏢，飛鏢扎進木板的任何一點的概率的相等的。這樣，我們投擲多次，例如100次，然後我們統計這100次中，扎入不規則圖形內部的次數，假設為k，那麼，我們就可以用 k/100 * 1 近似估計不規則圖形的面積，例如100次有32次擲入圖形內，我們就可以估計圖形的面積為0.32平方米。
以上這個過程，就是Monte-Carlo演算法直觀應用例子。
非形式化地說，Monte-Carlo演算法泛指一類演算法。在這些演算法中，要求解的問題是某隨機事件的概率或某隨機變數的期望。這時，通過「實驗」方法，用頻率代替概率或得到隨機變數的某些數字特徵，以此作為問題的解。
上述問題中，如果將「投擲一次飛鏢並擲入不規則圖形內部」作為事件，那麼圖形的面積在數學上等價於這個事件發生的概率（稍後證明），為了估計這個概率，我們用多次重複實驗的方法，得到事件發生的頻率 k/100 ，以此頻率估計概率，從而得到問題的解。
從上述可以看出，Monte-Carlo演算法區別於確定性演算法，它的解不一定是準確或正確的，其準確或正確性依賴於概率和統計，但在某些問題上，當重複實驗次數足夠大時，可以從很大概率上（這個概率是可以在數學上證明的，但依賴於具體問題）確保解的準確或正確性，所以，我們可以根據具體的概率分析，設定實驗的次數，從而將誤差或錯誤率降到一個可容忍的程度。
上述問題中，設總面積為S，不規則圖形面積為s，共投擲n次，其中擲在不規則圖形內部的次數為k。根據伯努利大數定理，當試驗次數增多時，k/n依概率收斂於事件的概率s/S。下面給出嚴格證明：
上述證明從數學上說明用頻率估計不規則圖形面積的合理性，進一步可以給出誤差分析，從而選擇合適的實驗次數n，以將誤差控制在可以容忍的範圍內，此處從略。
從上面的分析可以看出，Monte-Carlo演算法雖然不能保證解一定是準確和正確，但並不是「撞大運」，其正確性和準確性依賴概率論，有嚴格的數學基礎，並且通過數學分析手段對實驗加以控制，可以將誤差和錯誤率降至可容忍範圍。
3
Monte-Carlo演算法的數理基礎
這一節討論Monte-Carlo演算法的數理基礎。
首先給出三個定義：優勢，一致，偏真。這三個定義在後面會經常用到。
1) 設p為一個實數，且0.5
2) 如果對於同一實例，某Monte-Carlo演算法不會給出不同的解，則認為該演算法時一致的。
3) 如果某個解判定問題的Monte-Carlo演算法，當返回true時是一定正確的。則這個演算法時偏真的。注意，這裡沒有定義「偏假」，因為「偏假」和偏真是等價的。因為只要互換演算法返回的true和false，「偏假」就變成偏真了。
下面，我們討論Monte-Carlo演算法的可靠性和誤差分析。
總體來說，適用於Monte-Carlo演算法的問題，比較常見的有兩類。一類是問題的解等價於某事件概率，如上述求不規則圖形面積的問題；另一類是判定問題，即判定某個命題是否為真，如主元素存在性判定和素數測試問題。
先來分析第一類。對於這類問題，通常的方法是通過大量重複性實驗，用事件發生的頻率估計概率。之所以能這樣做的數學基礎，是伯努利大數法則：事件發生的頻率依概率收斂於事件的概率p。這個法則從數學生嚴格描述了頻率的穩定性，直觀意義就是當實驗次數很大時，頻率與概率偏差很大的概率非常小。此類問題的誤差分析比較繁雜，此處從略。有興趣的朋友可以參考相關資料。
接著，我們分析第二類問題。這裡，我們只關心一致且偏真的判定問題。下面給出這類問題的正確率分析：
由以上分析可以看到，對於一致偏真的Monte-Carlo演算法，即使調用一次得到正確解的概率非常小，通過多次調用，其正確率會迅速提高，得到的結果非常可靠。例如，對一個q為0.5的問題，假設p僅為0.01，通過調用1000次，其正確率約為0.9999784，幾乎可以認為是絕對準確的。重要的是，使用Monte-Carlo演算法解判定問題，其正確率不隨問題規模而改變，這就使得僅需要損失微乎其微的正確性，就可以將演算法複雜度降低一個數量級，在後面中可以看到具體的例子。
4
應用實例一
使用Monte-Carlo演算法計算定積分
計算定積分是金融、經濟、工程等領域實踐中經常遇到的問題。通常，計算定積分的經典方法是使用Newton-Leibniz公式：
這個公式雖然能方便計算出定積分的精確值，但是有一個局限就是要首先通過不定積分得到被積函數的原函數。有的時候，求原函數是非常困難的，而有的函數，如f(x) = (sinx)/x，已經被證明不存在初等原函數，這樣，就無法用Newton-Leibniz公式，只能另想辦法。
下面就以f(x) = (sinx)/x為例介紹使用Monte-Carlo演算法計算定積分的方法。首先需要聲明，f(x) = (sinx)/x在整個實數域是可積的，但不連續，在x = 0這一點沒有定義。但是，當x趨近於0其左右極限都是1。為了嚴格起見，我們補充定義當x = 0時f(x) = 1。另外為了需要，這裡不加證明地給出f(x)的一些性質：補充x = 0定義後，f(x)在負無窮到正無窮上連續、可積，並且有界，其界為1，即f(x)
下面開始介紹Monte-Carlo積分法。為了便於比較，在本節我們除了介紹使用Monte-Carlo方法計算定積分外，同時也探討和實現數值計算中常用的插值積分法，並通過實驗結果數據對兩者的效率和精確性進行比較。
1.插值積分法
我們知道，對於連續可積函數，定積分的直觀意義就是函數曲線與x軸圍成的圖形中，y>0的面積減掉y
圖2、梯形插值
如圖2所示，藍色部分是x1到x2積分的精確面積，而在梯形插值中，用橙色框所示的梯形面積近似估計積分值。
顯然，梯形法則的效果一般，而且某些情況下偏差很大，於是，有人提出了一種改進的方法：首先將積分區間分段，然後對每段計算梯形插值再加起來，這樣精度就大大提高了。並且分段越多，精度越高。這就是復化梯形法則。
除了梯形插值外，還有許多插值積分法，比較常見的有Sinpson法則，當然對應的也有復化Sinpson法則。下面給出四種插值積分的公式：
下面是四種插值積分法的程序代碼，用C#編寫。
2.Monte-Carlo積分法
我們知道，求定積分的直觀意義就是求面積，所以，用Monte-Carlo求積分的原理就是通過模擬統計方法求解面積。即通過向特定區域隨機產生大量點，然後統計點落在函數區域內的頻率，以此頻率估計面積，從而得到積分值。下面給出Monte-Carlo求取積分的演算法程序。
3.積分法的測試與比較
圖3、積分法測試結果
表1、積分方法實驗結果
首先看時間效率。當頻度較低時，各種方法沒有太多差別，但在1000萬級別上復化梯形和復化Sinpson相差不大，而Monte-Carlo演算法的效率快一倍。
而從準確率分析，當頻度較低時，幾種方法的誤差都很大，而隨著頻度提高，插值法要遠遠優於Monte-Carlo演算法，特別在1000萬級別時，Monte-Carlo法的相對誤差是插值法的的近萬倍。總體來說，在數值積分方面，Monte-Carlo方法效率高，但準確率不如插值法。
5
應用實例二
在O(n)複雜度內判定主元素
這次，我們看一個判定問題。問題是這樣的：在一個長度為n的數組中，如果有超過[n/2]的元素具有相同的值，那麼具有這個值的元素叫做數組的主元素。現在要求給出一種演算法，在O(n)時間內判定給定數組是否存在主元素。
如果採用確定性演算法，由於最壞情況下要搜索n/2次，而每次要比較的次數為O(n)量級，這樣，演算法的複雜度就是O(n^2)，不可能在O(n)時間內完成。所以我們只好換一種思路：不是要一個一定正確的結果，而只需要結果在很大概率上正確就行。我們可以這樣做：
圖4、Monte-Carlo法判定主元素
上述演算法，就是用Monte-Carlo思想求解主元素判定問題的過程。由於閾值N是一個給定的常數，不隨規模變化而變化，所以這個演算法的時間複雜度為O(n)，符合題設要求。但這個演算法給出的解並不是100%正確的，正確率和N有關。N設得過大，影響效率，N太小，正確率太低，那麼到底N設多大合適呢。這就要對演算法進行概率分析。
首先，這個演算法是一致且偏真的，證明很簡單，這裡從略。所以，如果數組中不存在主元素，則結果一定正確，而如果存在，調用一次得到正確結果的概率不低於1/2。由於偏真，在N次調用中只要返回一次True，就可以認為得到正確結果，所以，調用N此得到正確結果的概率不低於1 – (1/2)^N，可以看到，隨著N的增大，這個概率增加很快，只要調用10次，正確率就可以達到99.9%，重要的是，這個正確率和規模無關，即使數組的元素有1千萬億，只需調用10次，正確率依然是99.9%，這就體現出在數組很大時，Monte-Carlo方法的優勢。
下面是使用Monte-Carlo演算法進行主元素測試的C#程序示例。
程序很簡單，不做贅述。下面測試這個演算法。我們分別將閾值設為1、3、10，並且在每個閾值下測試100次，看看這個演算法的準確率如何。測試數組是[ 4, 5, 8, 1, 8, 4, 9, 2, 2, 2, 2, 2, 5, 7, 8, 2, 2, 2, 2, 2, 1, 0, 9, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 4, 7, 8, 2, 2, 2, 2, 2, 0, 1, 2, 2, 2, 2, 2 ]，其中存在主元素2。下面是測試結果：
圖5、Monte-Carlo演算法判定主元素實驗結果
測試數組有49個元素，主元素2有29個，比率為59%。從測試結果可以看出，即使閾值為1，正確率也高達84%，而僅僅為3的閾值就使正確率升高到98%，閾值為10時，100次測試全部正確。雖然理論上來說，閾值為10時有0.41^10=0.013%的概率給出錯誤判斷，但是筆者多次試驗，還沒有在閾值為10時得到錯誤結果。所以，Monte-Carlo方法求解判定問題，不論從理論上還是實踐中，都是不錯的方法。
另外一個與判定主元素類似的應用是素數判定問題，我們知道，對於尋找上百位的大素數，完全測試在時間效率上時不允許的。於是，結合費馬小定理使用Monte-Carlo法進行素數判定，是廣泛使用的方法。具體這裡不再詳述，感興趣的朋友可以參考相關資料。
6
應用實例三
分布未知的概率密度函數模擬
現在我們來看看Monte-Carlo演算法的第三種應用：模擬。在這種應用中，不再是用Monte-Carlo演算法求解問題，而是用來模擬難以解析描述的東西。問題是這樣的：
這個問題是實驗室一個師兄在開發Six Sigma軟體開發過程管理工具時遇到的一個實際需求，最終Y的概率密度函數將被用來計算分位點，從而進行過程式控制制。其中X可能是正態分布（高斯分布）、泊松分布、均勻分布或指數分布等。將多個不同分布的概率密度函數相加，得到的Y的分布式很難解析表示出來的，但如果是為了計算分位點，我們可以採取這樣一個策略：對於每一個X，產生若干符合其分布的點，帶入公式就得到若干符合Y分布的點，然後分段計算頻率，從而模擬出Y的分布，這些模擬點也可以用於分位點計算。這就是Monte-Carlo模擬的思想。
下面我們實現這個演算法，這裡的X我們僅給出最常用的正態分布，如果要實現其他分布，只要編寫相應的隨機點發生器就可以了。由於C#中只能產生符合均勻分布的隨機數，所以我們需要一種演算法，將均勻分布的隨機數轉為正態分布隨機數。這種演算法很多，Marc Brysbaert在1991年發表的《Algorithms for randomness in the behavioral sciences: A tutorial》一文中，共總結了5種將均勻分布隨機數轉為正態分布的隨機數的演算法，這裡筆者用到的是Knuth在1981年提出的一種演算法。這個演算法是將符合u(0,1)均勻分布的隨機點轉換為符合N(0,1)標準正態分布的隨機點p，由概率知識可知，要轉為符合N(e,v)的一般正態分布，只需進行p*v+e即可。下面是這個演算法：
下面是根據這個演算法，使用C#編寫的正態分布隨機點發生器：
接著是利用這個正態分布發生器獲得X的隨機值，並計算出Y的隨機值的代碼。也就是Y的隨機點發生器：
這樣，我們就可以產生任意多個符合Y分布的隨機點，從而藉此模擬Y的概率密度分布。
接著，我們測試一下這個模擬程序的效果，首先我們將初始值設為僅有一個符合標準正態分布的X，這樣Y=X，我們看看直接模擬一個標準正態分布的效果。這裡，我們產生100萬個隨機點。
圖6、使用Monte-Carlo演算法模擬標準正態分布
可以看到，模擬效果基本令人滿意。接下來，我們實際應用這個程序模擬一個分布未知的Y，其中Y = 15 + 2*N(2,8) + 5*N(-10,9) + 7*N(0,0.5)。模擬結果如下：
圖7、使用Monte-Carlo演算法模擬未知分布
有了符合Y分布的大量隨機點以及頻率統計，就可以隨心所欲繪出分布模擬圖，並進行分位點計算。這樣就用Monte-Carlo演算法解決了本節開頭提到的問題。
7
總結
本文首先通過一個不規則圖形面積計算的例子直觀介紹了Monte-Carlo演算法，然後給出了Monte-Carlo演算法在應用過程中需要了解的數理基礎。然後大篇幅介紹了三個應用：計算、判定和模擬。
總體來說，當需要求解的問題依賴概率時，Monte-Carlo方法是一個不錯的選擇。但這個演算法畢竟不是確定性演算法，在應用過程中需要冒一定「風險」，這就要求不能濫用這個演算法，在應用過程中，需要對其準確率或正確率進行數理分析，合理設計實驗，從而得到良好的結果，並將風險控制在可容忍的範圍內。
其實，不確定性演算法不只Monte-Carlo一種，Sherwood演算法、Las Vegas演算法和遺傳演算法等也是經典的不確定演算法。在很多問題上，不確定性演算法具有很好大的應用價值。有興趣的朋友可以參考相關資料。
參考文獻
[1] 孫海燕，周夢等 著，應用數理統計。北京航空航天大學出版社，2008.8
[2] 盛驟，謝式千，潘承毅 著，概率論與數理統計。高等教育出版社，2006.12
[3] David Kincaid，WardCheney 著，王國榮等 譯，數值分析(原書第三版)。機械工業出版社，2005.9
[4] Thomas H. Cormen等 著，演算法導論(第二版，影印版)。高等教育出版社，2002.5
[5] 王曉東 著，計算機演算法設計與分析。電子工業出版社，2001.1
[6] Marc Brysbaert，Algorithms for randomness in the behavioral sciences: A tutorial。Behavior Research Methods, Instruments, & Computers 1991, 23 (1) 45-60
[7] Patrick Smacchia 著，施凡等 譯，C#和.NET2.0 平台、語言與框架。2008.1
[8] Google。www.google.com
[9] Wikipedia。www.wikipedia.org
請您繼續閱讀更多來自 金融大數據研究 的精彩文章:TAG:金融大數據研究 |
※Material Design 控制項之 Toolbar 非完全解析
※Part3：About polite 解析
※Android 熱修復 Tinker Gradle Plugin解析
※Shoo-Career深度解析美國求職面試中的Case Interview
※UFC 213：副主賽 Yoel Romero vs.Robert Whittaker 賽前解析
※vuex所有核心概念完整解析State Getters Mutations Actions
※aws ec2的iam role深度解析
※MESA/Boogie Stowaway＆Highwire buffer單塊深度對比解析
※ReactiveCocoa源碼解析（一）Event與Observer代碼實現
※Unigine釋放顯卡性能測試Superposition，支持VR、8K解析度
※韓女團美妝亮點解析by.too cool for school
※vue之nextTick全面解析
※Flask 源碼解析：session
※《坦克世界》確認將為微軟Xbox Project Scorpio天蠍主機提供4K解析度
※崩壞學園2redcloak怎麼樣 redcloak實用性解析
※「乾貨」神經機器翻譯全流程解析，one-shot 和 zero-shot 學習成亮點
※JavaScript 中的 this 解析——兄弟連教育
※《W兩個世界》OST Part.4— 「Remember」MV及歌詞解析
※Coreldraw詳細解析女士T恤的款式圖畫法