機率算法中的Monte carlo算法

 最近在和同學讨論研究Six Sigma（六西格瑪）軟體開發方法及CMMI相關問題時，遇到了需要使用Monte-Carlo算法模拟分布未知的多元一次機率密度分布問題。于是花了幾天時間，通過查詢相關文獻資料，深入研究了一下Monte-Carlo算法，并以實際應用為背景進行了一些實驗。

      在研究和實驗過程中，發現Monte-Carlo算法是一個非常有用的算法，在許多實際問題中，都有用武之地。目前，這個算法已經在金融學、經濟學、工程學、實體學、計算科學及計算機科學等多個領域廣泛應用。而且這個算法本身并不複雜，隻要掌握機率論及數理統計的基本知識，就可以學會并加以應用。由于這種算法與傳統的确定性算法在解決問題的思路方面截然不同，作為計算機科學與技術相關人員以及程式員，掌握此算法，可以開闊思維，為解決問題增加一條新的思路。

      基于以上原因，我有了寫這篇文章的打算，一是回顧總結這幾天的研究和實驗，加深印象，二是和朋友們分享此算法以及我的一些經驗。

      這篇文章将首先從直覺的角度，介紹Monte-Carlo算法，然後介紹算法基本原理及數理基礎，最後将會和大家分享幾個基于Monte-Carlo方法的有意思的實驗。是以程式将使用C#實作。

      閱讀本文需要有一些機率論、數理統計、微積分和計算複雜性的基本知識，不過不用太擔心，我将盡量避免過多的數學描述，并在适當的地方對于用到的數學知識進行簡要的說明。

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算法雖然不能保證解一定是準确和正确，但并不是“撞大運”，其正确性和準确性依賴機率論，有嚴格的數學基礎，并且通過數學分析手段對實驗加以控制，可以将誤差和錯誤率降至可容忍範圍。

Monte-Carlo算法的數理基礎

      這一節讨論Monte-Carlo算法的數理基礎。

      首先給出三個定義：優勢，一緻，偏真。這三個定義在後面會經常用到。

      1) 設p為一個實數，且0.5<p<1。如果一個Monte-Carlo方法對問題任一執行個體的得到正确解的機率不小于p，則該算法是p正确的，且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算法解判定問題，其正确率不随問題規模而改變，這就使得僅需要損失微乎其微的正确性，就可以将算法複雜度降低一個數量級，在後面中可以看到具體的例子。

應用執行個體一：使用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)| <= 1，當且僅當x = 0時f(x) = 1。

      下面開始介紹Monte-Carlo積分法。為了便于比較，在本節我們除了介紹使用Monte-Carlo方法計算定積分外，同時也探讨和實作數值計算中常用的插值積分法，并通過實驗結果資料對兩者的效率和精确性進行比較。

1、插值積分法

      我們知道，對于連續可積函數，定積分的直覺意義就是函數曲線與x軸圍成的圖形中，y>0的面積減掉y<0的面積。那麼一種直覺的數值積分方法是通過插值方法，其中最簡單的是梯形法則：用以f(a)和f(b)為底，x軸和f(a)、f(b)連線為腰組成的梯形面積來近似估計積分。如下圖所示。

圖2、梯形插值

      如圖2所示，藍色部分是x1到x2積分的精确面積，而在梯形插值中，用橙色框所示的梯形面積近似估計積分值。

      顯然，梯形法則的效果一般，而且某些情況下偏差很大，于是，有人提出了一種改進的方法：首先将積分區間分段，然後對每段計算梯形插值再加起來，這樣精度就大大提高了。并且分段越多，精度越高。這就是複化梯形法則。

      除了梯形插值外，還有許多插值積分法，比較常見的有Sinpson法則，當然對應的也有複化Sinpson法則。下面給出四種插值積分的公式：

      下面是四種插值積分法的程式代碼，用C#編寫。

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82

using System; using System.Collections.Generic; using System.Linq; using System.Text; namespace MonteCarlo.Integration { /// <summary> /// 數值法求積分 /// 被積函數為 f(x) = (sin x)/x /// </summary> public class NumericalIntegrator { /// <summary> /// 梯形法則求積分 /// 積分公式為：((b - a) / 2) * [f(a) + f(b)] /// </summary> /// <param name="a">積分下限</param> /// <param name="b">積分上限</param> /// <returns>積分值</returns> public static double TrapezoidalIntegrate( double a,  double b) { return ((b - a) / 2) * (Math.Sin(a) / a + Math.Sin(b) / b); } /// <summary> /// 複化梯形法則求積分 /// 積分公式為：累加((xi - xi-1) / 2) * [f(xi) + f(xi-1)]  (i=1,2,...,n) /// </summary> /// <param name="a">積分下限</param> /// <param name="b">積分上限</param> /// <param name="n">分段數量</param> /// <returns>積分值</returns> public static double ComplexTrapezoidalIntegrate( double a,  double b,  int n) { double result = 0; for ( int i = 0; i < n; i++) { double xa = a + i * (b - a) / n; //區間積分下限 double xb = xa + (b - a) / n; //區間積分上限 result += ((xb - xa) / 2) * (Math.Sin(xa) / xa + Math.Sin(xb) / xb); } return result; } /// <summary> /// Sinpson法則求積分 /// 積分公式為：((b - a) / 6) * [f(a) + 4 * f((a + b) / 2) + f(b)] /// </summary> /// <param name="a">積分下限</param> /// <param name="b">積分上限</param> /// <returns>積分值</returns> public static double SinpsonIntegrate( double a,  double b) { return ((b - a) / 6) * (Math.Sin(a) / a + 4 * (Math.Sin(a + b) / (2 * (a + b))) + Math.Sin(b) / b); } /// <summary> /// 複化Sinpson法則求積分 /// 積分公式為：累加(h / 3) * [f(x2i-2) + 4*(f(x2i-1)) + f(x2i)]  (i=1,2,...,n/2 h = (b - a) / n) /// </summary> /// <param name="a">積分下限</param> /// <param name="b">積分上限</param> /// <param name="n">分段數量(必須為偶數)</param> /// <returns>積分值</returns> public static double ComplexSinpsonIntegrate( double a,  double b,  int n) { double result = 0; for ( int i = 0; i < n / 2 - 1; i++) { double xa = a + 2 * i * (b - a) / n; //區間積分下限 double xb = xa + (b - a) / n; //區間積分限中點 double xc = xb + (b - a) / n; //區間積分上限 result += ((b - a) / (3 * n) * (Math.Sin(xa) / xa + 4 * (Math.Sin(xb) / xb) + Math.Sin(xc) / xc)); } return result; } } }

2、Monte-Carlo積分法

      我們知道，求定積分的直覺意義就是求面積，是以，用Monte-Carlo求積分的原理就是通過模拟統計方法求解面積。即通過向特定區域随機産生大量點，然後統計點落在函數區域内的頻率，以此頻率估計面積，進而得到積分值。下面給出Monte-Carlo求取積分的算法程式。

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

using System; using System.Collections.Generic; using System.Linq; using System.Text; namespace MonteCarlo.Integration { /// <summary> /// Monte-Carlo法求積分 /// 被積函數為 f(x) = (sin x)/x /// </summary> public class MonteCarloIntegrator { /// <summary> /// 用Monte-Carlo法求解積分值 /// </summary> /// <param name="a">積分下限</param> /// <param name="b">積分上限</param> /// <param name="N">模拟次數</param> /// <returns>積分值</returns> public static double MonteCarloIntegrate( int a,  int b,  int N) { Random random =  new Random(); int positivePointCount = 0; //y >=0 區間内落入函數曲線内的點數目 int negativePointCount = 0; //y < 0區間内落入函數曲線内的點數目 //統計y >= 0區間點分布 for ( int i = 0; i < N; i++) { double xCoordinate = random.NextDouble(); //随機産生的x坐标 double yCoordinate = random.NextDouble(); //随機産生的y坐标 xCoordinate = a + (b - a) * xCoordinate; //将x規格化到相應積分區間 //yCoordinate = 1 * yCoordinate;//将y規格化到相應區間 if (Math.Sin(xCoordinate) / xCoordinate >= yCoordinate) { positivePointCount++; } } //統計y < 0區間點分布 for ( int i = 0; i < N; i++) { double xCoordinate = random.NextDouble(); //随機産生的x坐标 double yCoordinate = random.NextDouble(); //随機産生的y坐标 xCoordinate = a + (b - a) * xCoordinate; //将x規格化到相應積分區間 yCoordinate = -1 * yCoordinate; //将y規格化到相應區間 if (Math.Sin(xCoordinate) / xCoordinate <= yCoordinate) { negativePointCount++; } } double positiveFrequency = ( double )positivePointCount / ( double )N; //y >= 0區間内函數内點頻率 double negativeFrequency = ( double )negativePointCount / ( double )N; //y < 0區間内函數内點頻率 return (positiveFrequency - negativeFrequency) * ( double )(b - a); } } }

3、積分法的測試與比較

      下面對各種積分方法進行測試，對sinx/x在[1,2]區間上進行定積分。其中，我們分别對複化梯形和複化Sinpson法則做分段為10，10000，和10000000的積分測試。另外，對Monte-Carlo法也投點數也分為10，10000，和10000000。測試結果如下：

圖3、積分法測試結果

      為了分析偏差，我們必須給出一個精确值。但是現在我手頭沒有這個積分的精确值，不過1000萬次的梯形法則和Sinpson法則已經精确度很高了，是以這裡就以0.65932985作為基本，進行誤差分析。下面給出分析結果：

表1、積分方法實驗結果

      首先看時間效率。當頻度較低時，各種方法沒有太多差别，但在1000萬級别上複化梯形和複化Sinpson相差不大，而Monte-Carlo算法的效率快一倍。

      而從準确率分析，當頻度較低時，幾種方法的誤差都很大，而随着頻度提高，插值法要遠遠優于Monte-Carlo算法，特别在1000萬級别時，Monte-Carlo法的相對誤差是插值法的的近萬倍。總體來說，在數值積分方面，Monte-Carlo方法效率高，但準确率不如插值法。

應用執行個體二：在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 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

using System; using System.Collections.Generic; using System.Linq; using System.Text; namespace MonteCarlo.Detection { public class PrincipalElementDetector { /// <summary> /// 使用Monte-Carlo發探測主元素 /// </summary> /// <param name="elements">所有元素</param> /// <param name="N">門檻值</param> /// <returns>是否存在主元素</returns> public static bool DetectPrincipalElement(IList< int > elements, int N) { Random random =  new Random(); bool result =  false ; for ( int i = 0; i <= N; i++) { int index = random.Next(0, elements.Count - 1); int element = elements[index]; int count = 0; for ( int j = 0; j < elements.Count; j++) { if (element == elements[j]) { count++; } } if (count >= elements.Count / 2) { result =  true ; break ; } } return result; } } }

      程式很簡單，不做贅述。下面測試這個算法。我們分别将門檻值設為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法進行素數判定，是廣泛使用的方法。具體這裡不再詳述，感興趣的朋友可以參考相關資料。

應用執行個體三：分布未知的機率密度函數模拟

      現在我們來看看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#編寫的正态分布随機點發生器：

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

using System; using System.Collections.Generic; using System.Linq; using System.Text; namespace MonteCarlo.DistributingSimulation { public class NormalDistributingGenerator { /// <summary> /// 産生符合正态分布的随機數 /// 正态分布的期望為expectation,方差為variance /// </summary> /// <param name="expectation">期望</param> /// <param name="variance">方差</param> /// <param name="N">産生的數量</param> /// <returns>随機數序列</returns> public static IList< double > GenerateNDRNumber( double expectation,  double variance,  int N) { Random random =  new Random(); IList< double > randomList =  new List< double >(); for ( int i = 0; i < N; i++) { double u1, u2, v, z, a; do { u1 = random.NextDouble(); u2 = random.NextDouble(); v = 0.8578 * (2 * u2 - 1); z = v / u1; a = 0.25 * Math.Exp(2); if (a < 1 - u1) { break ; } }  while (a > 0.295 / u1 + 0.35 || a > -Math.Log(u1, Math.E)); randomList.Add(z * Math.Sqrt(variance) + expectation); } return randomList; } } }

      接着是利用這個正态分布發生器獲得X的随機值，并計算出Y的随機值的代碼。也就是Y的随機點發生器：

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

using System; using System.Collections.Generic; using System.Linq; using System.Text; namespace MonteCarlo.DistributingSimulation { public class DistributingSimulator { /// <summary> /// 模拟多個正态分布之和的分布情況，産生符合複合分布的随機點 /// y = a0 + a1*N(e1,v1) + ... + an*N(en,vn) /// N(e,v)表示期望為e，方差為v的正态分布 /// </summary> /// <param name="a">常數列</param> /// <param name="e">期望列</param> /// <param name="v">方差列</param> /// <param name="N">産生模拟點的個數</param> /// <returns>模拟點序列</returns> public static IList< double > Simulate(IList< double > a,IList< double > e,IList< double > v, int N) { IList< double > result =  new List< double >(); IList<IList< double >> randomLists =  new List<IList< double >>(); int count = a.Count - 1; //産生各個自變量的随機序列 for ( int i = 1; i <= count; i++) { randomLists.Add(NormalDistributingGenerator.GenerateNDRNumber(e[i], v[i], N)); } //帶入公式 for ( int j = 0; j < N; j++) { double y = 0; for ( int k = 1; k <= count; k++) { y += a[k] * randomLists[k - 1][j]; } y += a[0]; result.Add(y); } return result; } } }

      這樣，我們就可以産生任意多個符合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算法解決了本節開頭提到的問題。

總結

      本文首先通過一個不規則圖形面積計算的例子直覺介紹了Monte-Carlo算法，然後給出了Monte-Carlo算法在應用過程中需要了解的數理基礎。然後大篇幅介紹了三個應用：計算、判定和模拟。

      總體來說，當需要求解的問題依賴機率時，Monte-Carlo方法是一個不錯的選擇。但這個算法畢竟不是确定性算法，在應用過程中需要冒一定“風險”，這就要求不能濫用這個算法，在應用過程中，需要對其準确率或正确率進行數理分析，合理設計實驗，進而得到良好的結果，并将風險控制在可容忍的範圍内。

      其實，不确定性算法不隻Monte-Carlo一種，Sherwood算法、Las Vegas算法和遺傳算法等也是經典的不确定算法。在很多問題上，不确定性算法具有很好大的應用價值。有興趣的朋友可以參考相關資料。