簡明例析蒙特卡洛（Monte Carlo）抽樣方法

en.wikipedia.org

幾個基本分布

概率密度函數（PDF）描述了隨機變數在某個取值點取值的概率，累計分布函數（CDF）描述了隨機變數在 範圍內取值的概率，是概率密度函數的積分。

下圖給出了均勻分布（Uniform Distribution）的概率密度函數和累積分布函數（來自Uniform distribution (continuous)）。

均勻分布就是在所有取值點上取值概率都一樣的分布函數。在python中對均勻分布抽樣，可以先用下式取一個(0,1)之間的隨機數，然後映射到(a,b)之間。

x = numpy.random.rand()

如果取值區間可以簡單分為N段，每段的取值概率不同（有不同的權重），則可以先求出累積分布函數，將(0,1)之間的隨機數映射到累積分布函數上。累積分布函數上對應於該值的隨機變數的值就是符合概率分布的抽樣值。

均勻分布是一種最簡單的分布。另一種基本且應用非常廣泛的分布是正態分布。

正態分布（Normal Distribution），又稱高斯分布（Gauss Distribution），其概率分布函數和累積分布函數如下圖所示：（Normal distribution）

正態分布的概率密度函數解析式為：，其中 是正態分布的中心點， 是正態分布的方差。正態分布以中心點 左右對稱。

Naive Method

上述回顧了一些基本概念和常用分布函數，下面進入正題。

我們要解決的問題是：如何在不知道目標概率密度函數的情況下，抽取所需數量的樣本，使得這些樣本符合目標概率密度函數。這個問題簡稱為抽樣，是蒙特卡洛方法的基本操作步驟。

雖然我們不知道目標概率密度函數的具體形式，但是對任一確定的採樣點，我們可以衡量他的相對量。

題圖中 計算問題每個抽樣點可以表示為(x,y)，是二維的。為了方便，本文僅討論一維的採樣。討論過程中用到的一些思路有助於理解高維的情形。

為了類比二維 計算問題，我們構造一個一維問題：考慮一條一米長的繩子上被記號筆隨機塗了N段，我們想知道這N段一共有多長。

我們的解決方法是在（0,1）上做隨機採樣，可以很容易確定每一個樣本點是否被記號筆塗抹。那麼要知道被塗抹的長度總和，只需要數一數總共有多少個採樣點被記號筆塗抹，然後除以總樣本數。採樣點是否被記號筆塗抹的判斷如下所示：

def eval(x): if x > 0.21111 and x < 0.3232: return True if x > 0.4957 and x < 0.9191: return True return False

這個問題實在是太簡單了，所以我們把他稱為Naive Method.

下圖是採樣點的分布情況，綠色的採樣點是被記號筆塗抹過的。總採樣點數越多，最終計算的數值就越準確，這符合大數定理。

取捨演算法（Acceptance-Rejection Method）

上述Naive方法用均勻分布隨機得到的樣本都是有效樣本。這是因為上節所構造的記號筆塗抹問題本就適用於均勻分布。

對於一個一般性的問題，我們不能簡單的使用某一個已知分布直接得到所有的樣本。但是，我們可以先構造一個分布生成一堆隨機數，然後再在這些樣本中按照一定的方法挑選。這個構造的分布稱為建議分布函數（Proposed Distribution Function）。這種演算法被稱為取捨演算法（Acceptance-Rejection Method）。

取捨演算法不需要對概率密度函數歸一化。如果知道了所要衡量的函數值的範圍，就可以用一個覆蓋了目標概率密度函數的分布函數作為建議分布函數。

取捨演算法的步驟可以用如下代碼概括：

N=1000 #number of samples neededi = 1X = np.array([])while i < N: u = np.random.rand() x = (np.random.rand()-0.5)*8 res = u < eval(x)/ref(x) if res: X = np.hstack((X,x[res])) #accept ++i

我們的目標概率密度函數是eval(x)，建議分布函數是ref(x)。首先按照建議分布函數抽樣，然後取一個(0,1)之間的隨機數，判斷該隨機數是否小於eval(x)/ref(x)，即落在目標分布函數範圍內，如果滿足條件，該樣本被保留，否則進行下一輪測試。

建議分布函數可以是均勻分布，正態分布，也可以是其他分布。

下圖是抽樣結果的示意圖：（左邊使用的建議函數是均勻分布，右邊是正態分布）

重要性採樣（Importance Sampling）

重要性抽樣試圖通過使用和目標概率分布有相似形狀的建議分布函數以減小方差。但是如果選擇的概率密度函數不好，反而會增大方差。這裡不討論如何選擇較優的建議概率分布函數，只給出重要性採樣的基本原理和操作方法。

重要性採樣本質上還是取捨演算法。我們首先使用建議分布函數生成一堆樣本，然後計算每個樣本對目標概率分布函數的貢獻值，又稱為權重。權重越大，該樣本被選入的幾率就越高。直觀的看，如果建議概率分布中某個樣本的選擇概率較小，而目標分布函數中該樣本的選擇概率較大，則該樣本的權重就會比較大，說明建議概率分布中對該樣本的「重視程度」不夠。嚴重時，建議分布函數在某個區域的樣本值嚴重不足，會導致最終得到的樣本偏差較大，從而帶來更嚴重的方差。

重要性採樣的操作步驟如下所示：

N=100000M=5000x = (np.random.rand(N)-0.5)*16w_x = eval(x)/ref(x)w_x = w_x/sum(w_x)w_xc = np.cumsum(w_x) #accumulateX=np.array([])for i in range(M): u = np.random.rand() X = np.hstack((X,x[w_xc>u][0]))

其中，w_xc是對歸一化後的權重計算的累計分布概率。每次取最終樣本時，都會先隨機一個(0,1)之間的隨機數，並使用這個累計分布概率做選擇。樣本的權重越大，被選中的概率就越高。

下圖是分別使用均勻分布和正態分布函數作為建議分布函數所得的最終樣本分布圖。

梅特羅波利斯演算法（Metropolis Method）

上述方法都是基於取捨演算法，都必須事先生成較多採樣點備選。梅特羅波利斯演算法則通過構建一個穩態系統，採樣就是在狀態之間轉移的過程。該穩態系統基於馬爾可夫鏈，生成新樣本只依賴於上一個樣本和狀態，具有更好的適用性。

梅特羅波利斯演算法可以簡單描述為：

x = 0.1 #initialize x0 to be 0.1 for j in range(N[i]): u = np.random.rand() x_star = np.random.normal(x,10) A = min(1,eval(x_star)/eval(x)) if u < A: x = x_star X=np.hstack((X,x))

1）給定一個初始值作為當前值（示例中為0.1）

2）對於每個當前值x，根據建議分布函數隨機一個新樣本x_star

樣本接受率A = min(1, eval(x_star)/eval(x))

系統以一定的概率（rand() < A) 接受x_star成為新當前值

若未接受，當前樣本再次被接納，並再次成為當前值

梅特羅波利斯演算法要正常運行，必須使所構建的系統經過足夠多的狀態轉移後，收斂到一個平穩狀態。到達平穩狀態後，下一次轉移到某一個狀態（具體到這裡是採樣點）的概率就成為一個平穩概率分布（此時採樣點被選中的概率就接近於目標分布函數中該點的值）。而要構建這樣一個穩態系統，有一個充分非必要條件是（即細緻平衡條件）：

若正態分布作為建議分布函數 ，則滿足 . 由於 ，可以驗證 ，所以滿足細緻平衡條件，即最終收斂到穩態分布。

本示例中，梅特羅波利斯演算法都使用了正態分布作為建議分布函數。並分別使用了兩個概率密度函數作為目標概率分布做實驗，最終獲得的樣本數據分布如下：

目標概率分布一：

目標概率分布二：

其他

本文參考了 @秦春林THEGIBOOK 《全局光照技術》等相關材料

上述python示常式序和運行結果可前往如下鏈接下載：