蒙特卡洛算法及其實現

從今天開始要研究Sampling Methods，主要是MCMC算法。本文是開篇文章，先來了解蒙特卡洛算法。

 

 

Contents

 

   1. 蒙特卡洛介紹

   2. 蒙特卡洛的應用

   3. 蒙特卡洛積分

 

 

 

1. 蒙特卡洛介紹

 

   蒙特卡羅方法（Monte Carlo method），也稱統計模擬方法，是二十世紀四十年代中期由於科學技術的

   發展和電子計算機的發明，而被提出的一種以概率統計理論為指導的一類非常重要的數值計算方法。是指使

   用隨機數（或偽隨機數）來解決很多計算問題的方法。與它對應的是確定性算法。蒙特卡羅方法在金融工程

   學，宏觀經濟學，計算物理學（如粒子輸運計算、量子熱力學計算、空氣動力學計算）等領域應用廣泛。

 

   蒙特卡羅方法於20世紀40年代美國在第二次世界大戰中研制原子彈的“曼哈頓計划”計划的成員S.M.烏拉姆

   和J.馮·諾伊曼首先提出。數學家馮·諾伊曼用馳名世界的賭城—摩納哥的Monte Carlo—來命名這種方法，

   為它蒙上了一層神秘色彩。在這之前，蒙特卡羅方法就已經存在。1777年，法國數學家布豐提出用投針實驗

   的方法求圓周率，這被認為是蒙特卡羅方法的起源。

 

   另外，擬蒙特卡洛算法在近幾年也獲得迅速發展。這種方法是用確定性的超均勻分布代替蒙特卡洛算法中的

   隨機數序列，對於某些特定問題計算速度比普通的蒙特卡洛算法高幾百倍。

 

   由於產生隨機數的隨機性，當我們用N個隨機點以蒙特卡羅方法來求解具體的問題時，其計算得到近似解的誤

   差值有大有小，但是肯定有一個確定的平均值，即一些誤差大於此值，而其余誤差小於此值。鑒於此，顯然肯

   定存在這樣的N個點，使得誤差的絕對值不大於平均值。如果我們能夠構造這樣的點集，就可以對原有的方法

   進行較大的改進。擬蒙特卡羅方法就是至於此而提出的，它致力於構造其誤差比平均誤差顯著要好的那種點集，

   而其求解形式與蒙特卡羅方法一致，只不過所用的隨機數不一樣。用蒙特卡羅方法求解問題時，影響結果好壞

   的主要是隨機數序列的均勻性。而擬蒙特卡羅方法中的具有低偏差的一致分布點集較偽隨機數序列更為均勻，

   而且用擬蒙特卡羅方法求解得到的是真正的誤差，避免了蒙特卡羅方法得到概率誤差的缺陷。

   由此可見用擬蒙特卡羅方法求解問題的關鍵是如何找到一個均勻散布的點集。目前常用的點集有GLP點集(好格

   子點集，good lattice point set)、GP點集(好點集，good point set)、Halton點集及其變體、

   Hammersley點集等。

 

   蒙特卡洛方法的理論基礎是大數定律。大數定律是描述相當多次數重復試驗的結果的定律，根據這個定律知道

   樣本數量越多，其平均就越趨近於真實值。

 

 

 

2. 蒙特卡洛的應用

 

   最經典的應用就是利用蒙特卡洛算法求圓周率。代碼如下

 

代碼：

#include <bits/stdc++.h>

#define MAX_ITERS 1000000

using namespace std;

double Rand(double L, double R)
{
    return L + (R - L) * rand() * 1.0 / RAND_MAX;
}

double GetPi()
{
    srand(time(NULL));
    int cnt = 0;
    for(int i = 0; i < MAX_ITERS; i++)
    {
        double x = Rand(-1, 1);
        double y = Rand(-1, 1);
        if(x * x + y * y <= 1)
            cnt++;
    }
    return cnt * 4.0 / MAX_ITERS;
}

int main()
{
    for(int i = 0; i < 10; i++)
        cout << GetPi() << endl;
    return 0;
}

3. 蒙特卡洛積分

 

   關於蒙特卡洛求積分，可以先參照如下文章。

 

   鏈接：http://cos.name/2010/03/monte-carlo-method-to-compute-integration/ 

 

   接下來用蒙特卡洛積分求自然常數。這是2015年阿里的一道筆試題。

 

   首先考慮如下積分

 

      

 

   接下來分別用蒙特卡洛積分和牛頓萊布尼茲公式計算，在蒙特卡洛方法中樣本很多時，它們的值應該相等。

 

   利用蒙特卡洛方法，圖像大致如下

 

      

 

    上述積分的目的是求陰影部分的面積，所以先在所標矩形內取對隨機點，

    對於每一對，考察是否滿足如下條件

 

         

 

    假設滿足上述條件的點有個，而全部的點有個，所以得到近似公式為

 

         

 

    而依據牛頓萊布尼茲公式可以得到

 

         

 

    這兩種方法結果應該是相等的，即有

 

          

 

    接下來寫寫代碼吧！

 

代碼：

 

#include <bits/stdc++.h>

#define MAX_ITERS 10000000

using namespace std;

struct Point
{
    double x, y;
};

double Rand(double L, double R)  
{  
    return L + (R - L) * rand() * 1.0 / RAND_MAX;  
} 

Point getPoint()
{
    Point t;
    t.x = Rand(1.0, 2.0);
    t.y = Rand(0.0, 1.0);
    return t;
}

double getResult()
{
    int m = 0;
    int n = MAX_ITERS;
    srand(time(NULL));
    for(int i = 0; i < n; i++)
    {
        Point t = getPoint();
        double res = t.x * t.y;
        if(res <= 1.0)
            m++;
    }
    return pow(2.0, 1.0 * n / m);
}

int main()
{
    for(int i = 0; i < 20; i++)
        cout << fixed << setprecision(6) << getResult() << endl;
    return 0;
}

 

觀察一下運行結果，效果還是不錯的。如下圖