<style> .red { color: red; } .blue{ color:blue; } </style> # 蒙地卡羅方法緣由 <img src="https://hackmd.io/_uploads/ryLLhs-Yh.png" height =400; style="display:block; margin-left: 0;auto; width:100%" /> <figcaption>(戰爭促進了科學的發展，但卻也奪走了無數條生命...)</figcaption> <br> 1940年代，美國新墨西哥州內的洛斯阿拉莫斯國家實驗室( Los Alamos National Laboratory ) 的科學家們正秘密地進行和核武器相關的研究計畫。當時他們正在研究核武器核心中的中子擴散問題(Neutron Diffussion)，儘管他們已經掌握了多數必要的資料(例如中子在和原子核發生碰撞前平均移動距離等...... )，然而該問題其中涉及的計算量過於龐大，無法用 deterministic 的數學方法來解決，然而 John von Neumann, Stanisław Marcin Ulam 和 Nicholas Constantine Metropolis 等人卻讓這件事情有了進展...。 當時擅長研究隨機問題的 Stanisław Marcin，在某次聽聞了物理學家Nicholas Constantine Metropolis 等人的會議報告後，便和同時也在場的 John von Neumann 討論起了如何用"隨機模擬"的方法去解決該問題。 之後 John von Neumann 和 Nicholas Constantine Metropolis 等人更是使用史上第一台電子數位電腦 ENINC，配合該方法，克服了核分裂問題所涉及的龐大計算量，此方法之後更進一步用於製造氫彈的研究中，可以說是"曼哈頓計畫"中最核心的模擬方法之一!!! 但因為這是秘密研究，所以就像特務一樣，必須有個代號來稱呼 :sunglasses:，而 Ulam 的叔叔常常跟親戚借錢去摩納哥的蒙地卡羅賭場賭博(~~想必是常常輸錢，不然怎麼會跟親戚借錢~~)且這個方法又跟機率有關，於是 Nicholas Constantine Metropolis 惡趣味地將其取名為 "蒙地卡羅"，這便是"蒙地卡羅方法"的由來。 之後在核武器的計畫結束前，更是有段和這方法有關的小插曲，物理學家費米(Enrico Fermi)在蒙地卡羅方法發明後才參與核武器的相關計畫，然而他加入後表示在十幾年前他已使用類似的方法，所以費米若是在一開始或早點加入，應該會更早完成計算 XDD 時至今日，蒙地卡羅方法已不侷限於核分裂的研究中，該方法廣泛地用於金融工程，機器學習，計算物理等多個領域中，接著就讓我們來體驗一下這神奇的方法吧~~~ # 瘋狂投針的人... 在開始介紹蒙地卡羅方法之前，讓我們再提一個人，那就是十八世紀的法國博物學家布豐。 某天他不知道發生什麼事了，對著平行等距的木地板投了一根長度小於木板寬度的針，接著再去看這根針跟木地板組合出來的平行線有沒有相交，然後他重複了此行為很多次，並去計算針和平行線相交的機率。 <img src="https://hackmd.io/_uploads/HJJj2lXY3.png" height =200; style="display:block; margin:auto; width:60%" /> <figcaption>(假設投了 n 次針，有 h 次與線相交，則會發現當 n 足夠大時， 2Ln/th 趨近於 Pi)</figcaption> </br> 筆者我覺得這件事本身其實蠻瘋狂的，誰吃飽沒事會一直去投針= =，然後去計算針和木板組出的平行線相交的次數阿 !!! :thinking_face: 。重點是，隨著投著的次數越多次，算出來的機率值竟然和 $\pi$ 有關係，因而求得了 $\pi$ 的近似值，只能說天才和瘋子只有一線之隔 XD ![](https://hackmd.io/_uploads/BJgr0xQYn.png) # 大數法則 由此可見，大至核分裂，小至投針問題，都能看見機率的蹤跡。 實際上，這兩件毫不相干的事情背後都涉及到了下方的"大數法則"，也就是蒙地卡羅方法背後的原理: :::info $Theorem(Strong\ Law \ of\ Large \ Number)$ : &nbsp;&nbsp;&nbsp;&nbsp;Let $\{X_n\}$ be a sequence of $iid$ random variable with $E(|X_1|)<\infty$, then we have $$\frac{S_n}{n}\overset{a.s.}\to \ E(X_1)\ \ as \ \ n \to \infty, \ where \ S_n=\sum_{i=1}^{n}X_i.$$ ::: 我相信知道怎麼證明上述定理的人應該不多(~~例如Truncate手法，但不是數學系的學生不需要讀懂證明~~)，白話來說，就是隨著實驗的次數越多，計算出來的平均結果會越接近其期望值，最簡單的例子就是你投一枚公正的硬幣 N 次，當 N 夠大時，硬幣出現正面的次數除以 N 會靠近1/2， 不過拜於近代電腦快速發展所賜，我們現今已可以在個人電腦簡單快速地使用蒙地卡羅方法(~~不需要自己下去投針~~)，接下來就讓我們簡單地用 Python 加上此方法來估計 $\int_{0}^{1}x^ndx，n=1,2,3...$吧( 以 $x^2$ 為例 )。 # Toy Example: 估計 x^n 在 [0,1] 上的積分, n =1,2,3, ....(以 x^2 為例) ### Method1 撒豆子 <img src="https://hackmd.io/_uploads/Hk49rjSt3.png" height =300; style="display:block; margin:auto; width:60%" /> </br> 如上圖所示，令 $f(x)=x^2$. $x=1$ 和 $y=0$圍成的區域為 A。想像你跟布豐一樣，準備從一袋裝滿豆子的袋子拿出針來均勻地撒在 [0,1]X[0,1]，隨著撒著的豆子數量越多，就像你投擲一枚公正的硬幣一樣，你會發現 <font class='red'>$$落在區域A 的豆子數/所有的豆子數 \approx P(豆子落在區域A的機率)。$$</font> 而 P(豆子落在區域A的機率)可以很直覺的想成 $\frac{Area(A)}{Area([0,1]X[0,1])}= Area(A)$，(或是說$P(A)\equiv\mu(A)$ defines a probability measure on [0,1]x[0,1], where $\mu$ is Lebesgue measure and $A$ is a measurable subset in [0,1]x[0,1].) 因此我們投的豆子數量越多，越能得到 $Area(A)=\int_{0}^{1}x^2dx$ 的近似值. P.S. Let $\{X_i\}$ be a seq of iid r.v. with unif distribution on $[0,1] \times [0,1]$, $g(x) = 1_{A}(x), x \in[0,1]\times[0,1]$, indicator function of Area A, then by law of large number, $$\frac{\sum_{i=1}^{n}g(X_i)}{n}\overset{a.s.}\to \ E(g(X_1))=P(X_1\in A)=Area(A)\ \ as \ \ n \to \infty. $$ ```python= import random as rd import matplotlib.pyplot as plt import numpy as np # 取得落在 [0,1]x[0,1] 的點 Point (x,y) def getRandomPoint(): x = rd.random() y = rd.random() return (x,y) # 蒙地卡羅模擬，n:numberOfTrials, k: power of x^k, showPlot:布林值，決定是否要把圖畫出來 def Monte(n, k, showPlot = False): fig, ax = plt.subplots() #落在目標區域面積裡的次數 inTheArea = 0 #收集落在區域內的點，For 繪圖用 pointsToDraw = [] # 撒 n 次豆子 for i in range(n): # 豆子落地於[0,1]x[0,1] 的點 (x,y) = getRandomPoint() # 檢查是否有撒入目標區域面積內，有的話落在目標區域面積裡的次數++ if(y <= x**k ): inTheArea += 1 pointsToDraw.append((x,y)) print(f'x^{k} 於 [0,1] 區間面積估計為:', inTheArea/n,', numOfTrials:',n) if(showPlot): # 畫出 y= x^k 和 [0,1]X[0,1] 圍成的區域 x = np.linspace(0,1,1000) y = x**k ax.plot(x,y, marker='o', markersize=1, color = 'g') ax.plot(x, 0*x, marker='o', markersize=1, color = 'g') ax.plot(0*x+1, x, marker='o', markersize=1, color = 'g') # 將落在區域內的點畫出來 for point in pointsToDraw: ax.plot(point[0], point[1], marker='o', markersize=0.5, color = 'r') plt.show() return inTheArea/n # 使用蒙地卡羅方法估計 x^2 在 [0,1] 區間的積分 Monte(10000,2,False) ``` 若我們將畫圖的 flag 打開(改成 True)，則可以得到下面的圖，綠色的線為區域A的邊界，紅點代表落在區域A內的點。 ![](https://hackmd.io/_uploads/ryyfKGrK3.png) #### 執行結果 ![](https://hackmd.io/_uploads/BktgbSxFn.png) ### Method2 大數法則 (Law of Large Number) 令 $g(x) = x^2$, $X_i, i=1,2,3,...$ 為獨立同分布的變數且 $X_1\sim U(0,1),$ uniform distribution on [0,1]，則我們會有$E(|g(X_1)|) = E(g(X_1))<E(1)=1<\infty.$ 於是根據大數法則，我們知道 $$\frac{\sum_{i=1}^{n}g(X_i)}{n}\overset{a.s.}\to \ E(g(X_1))=\int_{0}^{1}x^2dx \ \ as \ \ n \to \infty. $$ ``` python import random as rd # 取得落在 [0,1] 區間的隨機點 x def get_Random_X(): return rd.random() # 蒙地卡羅模擬，n:numberOfTrials, k: power of x^k def Monte(n,k): sumOf_xWithPowerK = 0 for i in range(n): x = get_Random_X() sumOf_xWithPowerK += x**k print(f'x^{k} 於 [0,1] 區間面積估計為:', sumOf_xWithPowerK/n,', numOfTrials:',n) return sumOf_xWithPowerK/n # 執行蒙地卡羅方法 Monte(10000,2) ``` #### 執行結果 ![](https://hackmd.io/_uploads/rJu7boBY3.png) 可以發現上述兩種方法估計的東西些微不同，一個方法是去算落在區域內的個數，另外一個是去算有點類似黎曼合的Sum，但是都對 $\displaystyle\int_{0}^1x^2dx = 1/3$ 做了近似的估計。 # Observation... 1. 將方法 1 撒豆子落在區域 A 的判別方式 ($y<x^2$) 換成 ($x^2+y^2\leq1$)，則可以求得 $\displaystyle\frac{\pi}{4}$ 的近似值。同理，可以簡單地用類似的方式求橢圓$\displaystyle\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$的面積，三維單位球的體積等...。 2. 上述的範例用的隨機變數是 Unif distribution，實際上從大數法則定理的敘述可以發現蒙地卡羅方法選用的random variable 不一定要 Unif distribution， 3. 用蒙地卡羅計算高維度的積分概念和上述一維的計算是類似的，例如給定獨立同分布的Random Vector $\{X_n\}$，且 $X_1$ 的 pdf 是 $\phi(x)，g(x):\mathbb{R}^n \to\mathbb{R}$ 是 measurable 且 $E(g(X_1)) =\displaystyle\int_{\mathbb{R}^n}g(x)\phi(x)\ dx<\infty，$則根據大數法則，$$\frac{\sum_{i=1}^{n}g(X_i)}{n}\overset{a.s.}\to \ E(g(X_1))\ \ as \ \ n \to \infty. $$ 可以發現今天若是期望值的計算是若是非常複雜的，或是說 n 是個非常大的數字，例如氣體的分子數，則可以用蒙地卡羅方法去進行估計 :)) # 參考資料 1. https://foundation.nmns.edu.tw/writing/past.php?v=1&id=192 2. https://en.wikipedia.org/wiki/Monte_Carlo_method 3. https://sites.google.com/a/g2.nctu.edu.tw/unimath/2021-02/%E6%A9%9F%E7%8E%87%E9%97%9C%E6%96%BC%E5%B8%83%E8%B1%90%E6%8A%95%E9%87%9D-%E7%82%BA%E4%BB%80%E9%BA%BC%E5%8F%B2%E4%B8%8A%E5%AF%A6%E9%A9%97%E6%9C%80%E7%B2%BE%E6%BA%96%E7%9A%84%E6%95%B8%E5%AD%B8%E5%AE%B6%E6%8B%89%E6%92%92%E9%87%8C%E5%B0%BC%E8%A2%AB%E6%87%B7%E7%96%91%E4%BD%9C%E5%BC%8A%E5%91%A2 4. https://www.value-at-risk.net/generalizing-the-crude-monte-carlo-estimator/4