-
前言
MC方法的關鍵在於如何對想要的分布進行采樣，獲得采樣點。換句話說就是如何生成滿足指定分布的隨機數。在該系列文章中，我們有一個默認的假設就是已經有了一個能產生均勻分布隨機數的機制，不管它是硬件生成的真隨機數，還是算法模擬的偽隨機數。關於偽隨機數的生成算法，如線性同余法或者移位寄存器法，請參考文獻[2]相關章節。 -
MC方法基本介紹
在概率應用中，如果某個分布\(p(\mathbf{x})\)的解析表達式很復雜，導致涉及到該分布的計算(如計算某函數關於\(p\)的期望)難以得到解析結果，那么可以用\(m\)個關於 分布\(p(\mathbf{x})\)的采樣點 \(\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{m}\)來近似替代分布\(p(\mathbf{x})\)。比如在求期望\(E_{p(\mathbf{x})}[f(\mathbf{x})]\)時，可以用\(\sum_{i=1}^{m}f(\mathbf{x}_{i})\)近似替代，即：
\[E_{p(\mathbf{x})}[f(\mathbf{x})]\approx \sum_{i=1}^{m}f(\mathbf{x}_{i}) \]這種計算積分的方法就是MC。
另外，在計算積分\(\int_{V}f(\mathbf{x})d\mathbf(x)\)時，也可以利用MC方法得到近似解。
首先將積分式重寫為\(\int_{V}f(\mathbf{x})d\mathbf(x)=V\int_{V}f(\mathbf{x})\frac{1}{V}d\mathbf(x)=E_{p(\mathbf{x})}[f(\mathbf{x})]\)。
其中\(p(\mathbf{x})=\frac{1}{V}\)為在區域\(V\)中的均勻分布。所以計算積分\(\int_{V}f(\mathbf{x})d\mathbf(x)\)，首先在積分區域\(V\)中 均勻采樣 ，得到\(m\)個采樣點\(\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{m}\)，然后用求和\(V\sum_{i=1}^{m}f(\mathbf{x}_{i})/m\)代替原始積分，即
\[\int_{V}f(\mathbf{x})d\mathbf(x)\approx V\frac{\Sigma_{i=1}^{m}f(\mathbf{x}_{i})}{m} \] -
MC方法與固定網格法的優缺點
文獻
在一維情況下，假設采樣點數為\(m\)，固定網格法的誤差為\(m^{-1}\)，而MC方法的誤差為\(m^{-1/2}\)，固定網格法的精度更高。
在\(N\)維情況下，固定網格法想獲得\(m^{-1}\)的誤差，需要采樣點數為\(m^{N}\)。而MC方法的誤差永遠是\(m^{-1/2}\)，與問題維度無關。 -
鏈式模型的exact method
文獻 2.4節
鏈式模型概率分布的兩個特點：
- 可以用動態規划的方式求極值
- 可以用動態規划的方式計算邊緣分布，進而對模型進行采樣。
對高維概率分布的采樣可以通過依次求各個分量的邊緣分布並采樣得到。
比如對概率分布\(\pi(x_{1},x_{2},\cdots,x_{}{d})\)采樣，可以先求出邊緣分布\(\pi(x{d})=\int\pi(x_{1},x_{2},\cdots,x{d})dx_{1}d_{x2}\cdots,dx_{d-1}\)，然后對\(\pi(x{d})\)采樣，得到樣本點\(x_{d}^{0}\)。然后把\(x_{d}^{0}\)帶入到原分布\(\pi(x_{1},x_{2},\cdots,x_{d})\)，得到降了一維的分布\(\pi(x_{1},x_{2},\cdots,x_{d}^{0})\)。重復此過程直到對所有分量都進行采樣，便得到高維分布的一個采樣。如果分布\(\pi(x_{1},x_{2},\cdots,x_{}{d})\)剛好是鏈式結構的話，按照鏈式結構的反方向順序對各分量進行采樣可以充分利用動態規划算法計算各邊緣分布的歸一化常數(配分函數)，這種計算方法又叫 forward summation backward sampling。詳細的內容可以參見文獻 2.4.2節
-
Rejection Method
Rejection Method可以認為是從概率密度函數（Probability density function，pdf）的定義出發構造的一種比較高效的采樣方法。
假設\(p(\mathbf{x})\)為一pdf，由定義可知，\(\mathbf{x}\)落在\([\mathbf{x}_{0},\mathbf{x}_{0}+\Delta \mathbf{x}]\)區間的概率即為\(p(\mathbf{x})\)在該區間下覆蓋的面積。所以，如果我們能得到一個在\(p(\mathbf{x})\)覆蓋的面積（記做\(P\)）內的(二維)均勻分布的話，那么該分布的\(\mathbf{x}\)坐標的概率分布就滿足\(p(\mathbf{x})\)。
得到上述均勻分布的一個簡單方法是在一個包圍\(P\)的區域(記做\(Q\))內產生二維的均勻隨機樣本，然后只接受落在\(P\)內部的樣本，落在\(P\)外的樣本就丟棄（拒絕）。然后取這些隨機樣本的\(\mathbf{x}\)坐標，就可以得到\(p(\mathbf{x})\)對應的隨機抽樣點。
根據以上的討論可以發現，Rejection Method的效率與\(P\)和\(Q\)的面積之比有關，\(Q\)越大，則隨機樣本點落在\(P\)外部的概率也就越大，即被拒絕的概率也就越大。一個極端的例子是在整個二維平面內生成均勻分布隨機樣本點作為\(Q\)，此時，雖然理論上也能得到想要的抽樣點，但是任一采樣點被拒絕的概率都為1。
所以Rejection Method的關鍵點就在於選擇合適的\(Q\)，使其滿足以下要求：
- \(Q\)可以完全覆蓋\(P\)
- \(Q\)不會太大，以免影響效率
- \(Q\)內部的均勻分布容易產生
一般可以取\(Q\)同樣是一個函數\(q(\mathbf{x})=Mg(\mathbf{x})\)覆蓋的面積，其中\(g(\mathbf{x})\)是一個概率分布的pdf。此時，獲取\(Q\)內的二維均勻分布樣本點\((\mathbf{x},y)\)可以分兩步來進行：
- 得到\(\mathbf{x}\)分量的分布，根據前述討論，這是關於\(g(\mathbf{x})\)的分布乘以一個常量\(M\)。
- 得到\(y\)分量在\(Q\)內的均勻分布
其中第一步也要求\(g(\mathbf{x})\)的采樣可以較為容易的計算出來，比如通過cdf變換的方式。
綜上，利用Rejection Method對概率分布\(p(\mathbf{x})\)進行抽樣的流程可以歸納如下：
- 尋找一個概率分布\(g(\mathbf{x})\)，使其滿足\(Mg(\mathbf{x})\leq p(\mathbf{x})\)對所有\(\mathbf{x}\)成立(比如最簡單的均勻分布)
- 對\(g(\mathbf{x})\)進行采樣，得到一個采樣點\(\mathbf{x}_{0}\)
- 對[0,1]上的均勻分布進行采樣，得到采樣點\(y_{0}\)
- 如果\(y_{0}\leq\frac{p(\mathbf{x})}{Mg(\mathbf{x})}\)，則接受\(\mathbf{x}_{0}\)，否則拒絕並轉向步驟2。
其中步驟2&3相當於是得到\(Q\)的均勻分布采樣點，步驟4是判斷該采樣點是否落在\(P\)內
關於Rejection Method的數學證明，可以參考文獻，非常簡單。
示例代碼如下：
%matplotlib inline import numpy as np import matplotlib.pyplot as plt import seaborn as sns # use uniform distrubituion as proposal distrubituion m = 1000000 bandwidth = 10 x = np.random.rand(m) * bandwidth x.sort() # define the target distribution mu1 = bandwidth * 0.25 mu2 = bandwidth * 0.75 sigma1 = 0.5 sigma2 = 1 p = 0.4 / (np.sqrt(2 * np.pi) * sigma1**2) * np.exp(-0.5 * ( x - mu1)**2 / sigma1**2) + 0.5 / (np.sqrt(2 * np.pi) * sigma2**2) * np.exp( -0.6 * (x - mu2)**2 / sigma2**2) # compute the reject propability M = 10 r = p / M * bandwidth # generate the samples x2 = np.random.rand(m) accepted_samples_x = x[x2 < r] # plot n, bins, patches = plt.hist( accepted_samples_x, bins=100, normed=True, color=sns.desaturate("indianred", .8), label="sample distribution") plt.plot(x, p / max(p) * max(n), 'blue', label="target distribution",linewidth=4) plt.legend(fontsize=15) plt.show() -
Importance Sampling
importance sampling的目的是為了計算積分\(E_{\pi}[h(\mathbf{x})]=\int h(\mathbf{x})\pi(\mathbf{x})\)。和rejection method不同，importance sampling本身並不能也不期望得到分布\(\pi(\mathbf{x})\)的采樣，而是轉而去對一個很容易采樣的分布\(g(\mathbf{x})\)進行采樣，然后在計算積分的時候加考慮采樣點各自的權重，以逼近真實的積分。這樣做的好處是除了可以避免對復雜分布\(\pi(\mathbf{x})\)采樣以外，還能通過設計\(g(\mathbf{x})\)的形式，使得采樣點更多的出現在對積分貢獻大的區域，以提高計算的精度和效率。
基本的importance sampling原理可以表示如下：
\[E_{\pi}[h(\mathbf{x})]=\int h(\mathbf{x})\pi(\mathbf{x})=\int h(\mathbf{x})g(\mathbf{x})\frac{\pi(\mathbf{x})}{g(\mathbf{x})}=E_{g}[\omega h(\mathbf{x})] \]其中\(\omega=\frac{\pi(\mathbf{x})}{g(\mathbf{x})}\)為權重，\(g(\mathbf{x})\)是用來近似\(h(\mathbf{x})\pi(\mathbf{x})\)的分布。
-
結合Rejection Method和Importance Sampling
利用Rejection Method得到采樣分布並計算積分的缺點是要選擇合適的\(M\)才能得到較高的效率，因為平均每\(M\)個采樣點才會有一個被接受。但是\(M\)又不能選擇的太小以避免在某些區域出現\(Mg(\mathbf{x})< p(\mathbf{x})\)的情況。
而利用Importance Sampling計算積分時，雖然對測試分布沒有什么要求（這點和Rejection Method不太一樣，Rejection Method要求測試分布\(g(\mathbf{x})\)一定要滿足\(Mg(\mathbf{x})\leq p(\mathbf{x})\)），但是如果測試分布與目標分布的差別非常大，那么在計算權重時就會出現大多數采樣點的權重都非常小的情況，從而也會影響計算效率
為了解決以上問題，發展出了一種綜合Rejection Method和Importance Sampling的方法，即Rejection Control。
Rejection Control可以認為是結合了Rejection Method的Importance Sampling方法。假設需要對分布\(\pi(\mathbf{x})\)進行采樣，得到\(m\)個采樣點\(\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{m}\)，並用它計算某個積分\(E_{\pi}[h(\mathbf{x})]=\int h(\mathbf{x})\pi(\mathbf{x})\)。根據Importance Sampling的思路，可以選取一個測試分布\(g(\mathbf{x})\)，對其進行采樣，然后給每個采樣點\(\mathbf{x}_{j}\)賦一個權重\(w_{i}=\frac{\pi(\mathbf{x})}{g(\mathbf{x})}\)用於計算積分\(E_{\pi}[h(\mathbf{x})]\)。
與普通Importance Sampling不同的是，在采樣的時候，增加以下步驟：
- 以概率\(r_{i}=\textrm{min}\{1,\frac{w_{i}}{c}\}\)接受該采樣點，其中\(c\)是一固定常數
- 如果某個采樣點\(\mathbf{x}_{i}\)被接受，那么其權重\(w_{i}\)更新為\(\bar{w}_{i}=q_{c}\frac{w_{i}}{r_{i}}\)，其中\(q_{c}=\int \textrm{min}\{1,\frac{w}{c}\}g(\mathbf{x})d\mathbf{x}\)為對所有采樣值都相同的常數。如果只是計算積分的話，該常數可以忽略，因為根據中(2.4)式，對所有采樣點都相同的常量會被約掉。
Rejection Control方法也可以認為是Rejection Method的改良版，利用采樣點權重\(w\)來衡量測試分布與目標分布的差距，並據此更新測試分布。根據\(w_{i}=\frac{\pi(\mathbf{x})}{g(\mathbf{x})}\)可得
\[\bar{g}(\mathbf{x})=\frac{\pi(\mathbf{x})}{\bar{w}}=\frac{\pi(\mathbf{x})}{q_{c}\frac{w}{r}}=\frac{\pi(\mathbf{x})r}{q_{c}w}=\frac{\textrm{min}\{\pi(\mathbf{x}),\frac{\pi(\mathbf{x})w}{C}\}}{q_{c}w}=q_{c}^{-1}\textrm{min}\{g(\mathbf{x}),\frac{\pi(\mathbf{x})}{c}\} \]可見，當\(q_{c}=c=1\)時，\(g(\mathbf{x})\)在某些地方會等於\(\pi(\mathbf{x})\)。
參考文獻
[1] Monte Carlo Strategies in Scientific Computing
[2] Numerical recipes
[3] Monte Carlo Statistical Methods
[4] Machine learning: a probabilistic perspective