蒙特卡洛算法及其实现
從今天開始要研究Sampling Methods，主要是MCMC算法。本文是開篇文章，先來了解蒙特卡洛算法。
?
?
Contents
?
?? 1. 蒙特卡洛介紹
?? 2. 蒙特卡洛的應(yīng)用
?? 3. 蒙特卡洛積分
?
?
?
1. 蒙特卡洛介紹
?
?? 蒙特卡羅方法（Monte Carlo method），也稱統(tǒng)計(jì)模擬方法，是二十世紀(jì)四十年代中期由于科學(xué)技術(shù)的
?? 發(fā)展和電子計(jì)算機(jī)的發(fā)明，而被提出的一種以概率統(tǒng)計(jì)理論為指導(dǎo)的一類非常重要的數(shù)值計(jì)算方法。是指使
?? 用隨機(jī)數(shù)（或偽隨機(jī)數(shù)）來解決很多計(jì)算問題的方法。與它對應(yīng)的是確定性算法。蒙特卡羅方法在金融工程
?? 學(xué)，宏觀經(jīng)濟(jì)學(xué)，計(jì)算物理學(xué)（如粒子輸運(yùn)計(jì)算、量子熱力學(xué)計(jì)算、空氣動(dòng)力學(xué)計(jì)算）等領(lǐng)域應(yīng)用廣泛。
?
?? 蒙特卡羅方法于20世紀(jì)40年代美國在第二次世界大戰(zhàn)中研制原子彈的“曼哈頓計(jì)劃”計(jì)劃的成員S.M.烏拉姆
???和J.馮·諾伊曼首先提出。數(shù)學(xué)家馮·諾伊曼用馳名世界的賭城—摩納哥的Monte Carlo—來命名這種方法，
?? 為它蒙上了一層神秘色彩。在這之前，蒙特卡羅方法就已經(jīng)存在。1777年，法國數(shù)學(xué)家布豐提出用投針實(shí)驗(yàn)
???的方法求圓周率，這被認(rèn)為是蒙特卡羅方法的起源。
?
?? 另外，擬蒙特卡洛算法在近幾年也獲得迅速發(fā)展。這種方法是用確定性的超均勻分布代替蒙特卡洛算法中的
?? 隨機(jī)數(shù)序列，對于某些特定問題計(jì)算速度比普通的蒙特卡洛算法高幾百倍。
?
?? 由于產(chǎn)生隨機(jī)數(shù)的隨機(jī)性，當(dāng)我們用N個(gè)隨機(jī)點(diǎn)以蒙特卡羅方法來求解具體的問題時(shí)，其計(jì)算得到近似解的誤
?? 差值有大有小，但是肯定有一個(gè)確定的平均值，即一些誤差大于此值，而其余誤差小于此值。鑒于此，顯然肯
?? 定存在這樣的N個(gè)點(diǎn)，使得誤差的絕對值不大于平均值。如果我們能夠構(gòu)造這樣的點(diǎn)集，就可以對原有的方法
?? 進(jìn)行較大的改進(jìn)。擬蒙特卡羅方法就是至于此而提出的，它致力于構(gòu)造其誤差比平均誤差顯著要好的那種點(diǎn)集，
???而其求解形式與蒙特卡羅方法一致，只不過所用的隨機(jī)數(shù)不一樣。用蒙特卡羅方法求解問題時(shí)，影響結(jié)果好壞
?? 的主要是隨機(jī)數(shù)序列的均勻性。而擬蒙特卡羅方法中的具有低偏差的一致分布點(diǎn)集較偽隨機(jī)數(shù)序列更為均勻，
?? 而且用擬蒙特卡羅方法求解得到的是真正的誤差，避免了蒙特卡羅方法得到概率誤差的缺陷。
?? 由此可見用擬蒙特卡羅方法求解問題的關(guān)鍵是如何找到一個(gè)均勻散布的點(diǎn)集。目前常用的點(diǎn)集有GLP點(diǎn)集(好格
?? 子點(diǎn)集，good lattice point set)、GP點(diǎn)集(好點(diǎn)集，good point set)、Halton點(diǎn)集及其變體、
?? Hammersley點(diǎn)集等。
?
?? 蒙特卡洛方法的理論基礎(chǔ)是大數(shù)定律。大數(shù)定律是描述相當(dāng)多次數(shù)重復(fù)試驗(yàn)的結(jié)果的定律，根據(jù)這個(gè)定律知道
?? 樣本數(shù)量越多，其平均就越趨近于真實(shí)值。
?
?
?
2. 蒙特卡洛的應(yīng)用
?
?? 最經(jīng)典的應(yīng)用就是利用蒙特卡洛算法求圓周率。代碼如下
?
代碼：
1 #include <bits/stdc++.h> 2 3 #define MAX_ITERS 1000000 4 5 using namespace std; 6 7 double Rand(double L, double R) 8 { 9 return L + (R - L) * rand() * 1.0 / RAND_MAX; 10 } 11 12 double GetPi() 13 { 14 srand(time(NULL)); 15 int cnt = 0; 16 for(int i = 0; i < MAX_ITERS; i++) 17 { 18 double x = Rand(-1, 1); 19 double y = Rand(-1, 1); 20 if(x * x + y * y <= 1) 21 cnt++; 22 } 23 return cnt * 4.0 / MAX_ITERS; 24 } 25 26 int main() 27 { 28 for(int i = 0; i < 10; i++) 29 cout << GetPi() << endl; 30 return 0; 31 }3. 蒙特卡洛積分
?
?? 關(guān)于蒙特卡洛求積分，可以先參照如下文章。
?
?? 鏈接：http://cos.name/2010/03/monte-carlo-method-to-compute-integration/?
?
?? 接下來用蒙特卡洛積分求自然常數(shù)。這是2015年阿里的一道筆試題。
?
?? 首先考慮如下積分
?
??????
?
???接下來分別用蒙特卡洛積分和牛頓萊布尼茲公式計(jì)算，在蒙特卡洛方法中樣本很多時(shí)，它們的值應(yīng)該相等。
?
?? 利用蒙特卡洛方法，圖像大致如下
?
??????
?
??? 上述積分的目的是求陰影部分的面積，所以先在所標(biāo)矩形內(nèi)取對隨機(jī)點(diǎn)，
??? 對于每一對，考察是否滿足如下條件
?
?????????
?
??? 假設(shè)滿足上述條件的點(diǎn)有個(gè)，而全部的點(diǎn)有個(gè)，所以得到近似公式為
?
?????????
?
????而依據(jù)牛頓萊布尼茲公式可以得到
?
?????????
?
??? 這兩種方法結(jié)果應(yīng)該是相等的，即有
?
??????????
?
????接下來寫寫代碼吧！
?
代碼：
?
1 #include <bits/stdc++.h> 2 3 #define MAX_ITERS 10000000 4 5 using namespace std; 6 7 struct Point 8 { 9 double x, y; 10 }; 11 12 double Rand(double L, double R) 13 { 14 return L + (R - L) * rand() * 1.0 / RAND_MAX; 15 } 16 17 Point getPoint() 18 { 19 Point t; 20 t.x = Rand(1.0, 2.0); 21 t.y = Rand(0.0, 1.0); 22 return t; 23 } 24 25 double getResult() 26 { 27 int m = 0; 28 int n = MAX_ITERS; 29 srand(time(NULL)); 30 for(int i = 0; i < n; i++) 31 { 32 Point t = getPoint(); 33 double res = t.x * t.y; 34 if(res <= 1.0) 35 m++; 36 } 37 return pow(2.0, 1.0 * n / m); 38 } 39 40 int main() 41 { 42 for(int i = 0; i < 20; i++) 43 cout << fixed << setprecision(6) << getResult() << endl; 44 return 0; 45 }?
觀察一下運(yùn)行結(jié)果，效果還是不錯(cuò)的。如下圖
?
?????????????
?
總結(jié)
以上是生活随笔為你收集整理的蒙特卡洛算法及其实现的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 网络协议基础
- 下一篇： 在Eclipse中使用JUnit4进行单