蒙特卡洛方法

常見使用場景

機器學習中經常會遇到對復雜的分布做加和或積分，例如在貝葉斯方法中，往往要對參數做積分，\(P(t|X)=\int p(t|\theta)p(\theta|X)d\theta\)，頻率派中EM算法的E步也是一個求期望的過程，\(Q(\theta,\theta_{old})=\int p(Z|X,\theta_{old})log\,p(Z,X|\theta)dZ\)，這些積分或者期望往往都是intractable的。對於這個問題，我們可以使用變分來解決，變分是通過最優化方法尋找一個和\(p(\theta|X)\)或\(p(Z|X,\theta_{old})\)近似且好積分的分布。而在這篇博文中，我們將介紹sampling的方法，舉個簡單的例子：比如我們遇到這種形式\(E[f]=\int f(z)p(z)\,dz\)，如果我們從p(z)中sampling一個數據集\(z^{(l)}\)，然后再求個平均\(\hat f=\frac{1}{L}\sum_{l=1}^L f(z^{(l)})\)來近似\(f(z)\)的期望，so問題就解決了，關鍵是如何從\(p(z)\)中做無偏的sampling。現在的問題轉換為給定一個概率分布，我們如何在計算機中生成它的樣本。一般而言均勻分布的樣本是相對容易生成的。 通過線性同余發生器可以生成偽隨機數，這些偽隨機數序列的各種統計指標和均勻分布的理論計算結果非常接近。這樣的偽隨機序列就有比較好的統計性質，可以被當成真實的隨機數使用。而我們常見的概率分布，無論是連續的還是離散的分布，都可以基於均勻分布的樣本生成。

蒙特卡洛方法

1.Inverse transform sampling

我們知道，計算機本身是無法產生真正的隨機數的，但是可以根據一定的算法產生偽隨機數。最古老最簡單的莫過於 線性同余發生器：$$x_{n+1}=(ax_{n}+c);mod;m$$式子中的\(a\)和\(c\)是一些數學知識推導出的合適的常數。但是我們看到，這種算法產生的下一個隨機數完全依賴現在的隨機數的大小，而且當你的隨機數序列足夠大的時候，隨機數將出現重復子序列的情況。當然，理論發展到今天，有很多更加先進的隨機數產生算法出現，比如python數值運算庫numpy用的是Mersenne Twister等。但是不管算法如何發展，這些都不是本質上的隨機數，用馮諾依曼的一句話說就是：

Anyone who considers arithmetic methods of producing random digits is, of course, in a state of sin.

OK，根據上面的算法現在我們有了均勻分布的隨機數，但是如何產生滿足其他分布（比如高斯分布）下的隨機數呢？一種可選的簡單的方法是 Inverse transform sampling。它的原理是利用累積分布函數（CDF，cumulative distribution function）來處理。

• Inverse transform sampling
• 假設：已經一個生成器可以產生(0, 1)區間上的均勻分布隨機數
• 問題：從分布p(z)中采樣
• 算法：假設有分布\(p(x)\)的累計函數是\(P(x)\)，若\(y\)是(0, 1)上的均勻分布隨機變量，那么\(P^{-1}(y)\)的就是服從分布\(p(x)\)的隨機變量
• 缺陷：\(P^{-1}(y)\)不容易計算

舉個列子，我們想對高斯分布采樣，首先在 y 軸上產生（0,1）之間的均勻分布的隨機數，水平向右投影到高斯累計分布函數上，然后垂直向下投影到 x 軸，得到的就是高斯分布的一個實例。可見高斯分布的隨機數實際就是均勻分布隨機數在高斯分布的 CDF 函數下的逆映射。當然，在實際操作中，更有效的計算方法有Box-Muller_transform ，Ziggurat algorithm 等，這些方法 tricky and faster，沒有深入了解，這里也不多說了。

2.Rejection Sampling

Rejection Sampling和Importance Sampling都是基於proposal distribution的sampling，這兩種方法都是假設直接從分布p(z)采樣很困難，但可以找到一個相對更容易采用的分布q(z)，即proposal distribution。
假設分布\(p(z)=\frac{1}{Z}\hat p(z)\)，其中Z是\(p(z)\)中與z無關的一個因子。之所以要分離開來寫，是因為有時候\(p(z)\)中與z無關的部分可能比較復雜，難以計算。例如\(p(z|x)=\frac{p(x,z)}{p(x)},p(x)=\int p(x,z)dz\)，其中，\(p(x)\)與z無關且難以計算。因此，在大多數的sampling中，只需要利用與z有關部分即可進行采樣的算法，這樣就避開了對復雜的Z的計算。

• Rejection Sampling
• 假設：直接從分布p(z)中采樣困難
• 問題：從分布p(z)中采樣
• 算法：首先為p(z)找一個proposal distribution q(z)，而且q(z)必須是很容易進行sampling的。然后找到一個盡可能小的常數k，使得\(kq(z)\geq p(z)\)，對任意z成立。然后從q(z)中sample出一個數\(z_0\)；然后從均勻分布[0,\(kq(z_0)\)]中sample出另一個數\(u_0\)；這時候，平面上的點（\(z_0\)，\(u_0\)）是\(kq(z)\)下方區域中的均勻分布。如果\(u_0>\hat p(z_0)\)，則拒絕樣本\(z_0\)並重復前面步驟，否則接收\(z_0\)為符合分布p(z)的點。
• 缺陷：上述過程看出了k為什么要盡可能小。K越小，才能使\(z_0\)被拒絕的概率盡可能小，從而提高rejection sampling的效率，因此我們需要選擇一個合適的k 。維數越高，拒絕率越高，采樣效率越低。例如高維的球，可計算其測度主要集中在球的表面；而rejection sampling中，\(u_0>\hat p(z_0)\)的部分正是高維幾何體的表層。這就導致在高維情況下，有很高的拒絕率。

z從q(z)分布抽樣，而接受率是\(\frac{\hat p(z)}{kq(z)}\)，所以

\begin{split}
p(accept)&=\int {\hat p(z)/kq(z)}q(z)dz \
&= \frac{1}{k} \int \hat p(z)dz
\end{split}

3.Importance Sampling

Importance Sampling和Rejection Sampling類似都是假設對p(z)采樣比較困難，不過對於一個給定的z，卻可容易的計算其概率值p(z)。但是與Rejection Sampling不同的是，Importance Sampling不是求p(z)樣例，而是直接計算函數f(z)在該分布下的期望。因為p(z)本身采樣困難，所以我們還是得像Rejection sampling那樣，找到一個更容易采樣的分布q(z)，並且假設從q(z)采樣了L個樣本。那么：$$E[f(z)]=\int f(z)p(z)dz=\int \frac{f(z)p(z)}{q(z)}q(z)dz\approx \frac{1}{L}\sum_{l=1}^L\frac{p(z^l)}{q(z^l)}f(z^l) $$
其中\(r_l=\frac{p(z^l)}{q(z^l)}\)被成為importance weights。
再進一步假設：對分布p(z)的認識集中在p(z)的與z相關的部分\(\hat p(z)\)，其normalization constant \(Z_p\)還未知。同時也從q(z)中分離出一個常數\(Z_q\)，那么：

\begin{split}
E[f(z)]&=\int f(z)p(z)dz=\frac{Z_q}{Z_p}\int \frac{f(z)\hat p(z)}{\hat q(z)}q(z)dz\
&\approx \frac{Z_q}{Z_p}\frac{1}{L}\sum_{l=1}^L\frac{\hat p(z^l)}{\hat q(z^l)}f(z^l)
\end{split}

其中\(\hat r_l=\frac{\hat p(z^l)}{\hat q(z^l)}\)，同樣地，可以計算：

\[\frac{Z_p}{Z_q}=\frac{1}{Z_q}\int \hat p(z)dz=\int \frac{\hat p(z)}{\hat q(z)}q(z)dz\approx \frac{1}{L} \sum_{l=1}^L\hat r_l \]

於是最終得到：

\[E[f(z)]=\sum_{l=1}^Lw_lf(z^l) \]

其中，\(z^l\)是分布從q(z)采樣的L個樣本，而

\[w_l=\frac{\hat r^l}{\sum_{m=1}^L\hat r^m} \]

注意，在\(w_l\)的計算中，已經只需要分布p(z)的與z有關部分\(\hat p(z)\)。從而達到了避開Z的目的。

• Importance Sampling
• 假設：直接從分布p(z)中采樣困難
• 問題：計算函數f(z)在p(z)分布下的期望
• 算法：首先為p(z)找一個proposal distribution q(z)，而且q(z)必須是很容易進行sampling的並且和目標分布相似。然后從q(z)中sample出L個數\(z_l\)，帶入公式\(E[f(z)]=\sum_{l=1}^Lw_lf(z^l)\)求的f(z)在p(z)分布下的期望
• 缺陷：但是可惜的是，在高維空間里找到一個這樣合適的q非常難。即使有 Adaptive importance sampling 和 Sampling-Importance-Resampling(SIR) 的出現，要找到一個同時滿足容易抽樣並且和目標分布相似的proposal distribution，通常是不可能的！

馬爾可夫蒙特卡洛

由於Rejection sampling和Importance sampling這兩種方法在高維下都會失效。我們的目標還是對於給定的概率分布p(x),我們希望能有便捷的方式生成它對應的樣本。由於馬氏鏈能收斂到平穩分布，於是一個很的漂亮想法是：如果我們能構造一個轉移矩陣為P的馬氏鏈，使得該馬氏鏈的平穩分布恰好是p(x), 那么我們從任何一個初始狀態x0出發沿着馬氏鏈轉移, 得到一個轉移序列 x0,x1,x2,⋯xn,xn+1⋯,， 如果馬氏鏈在第n步已經收斂了，於是我們就得到了樣本xn,xn+1⋯。

這個絕妙的想法在1953年被Metropolis想到了，為了研究粒子系統的平穩性質，Metropolis考慮了物理學中常見的波爾茲曼分布的采樣問題，首次提出了基於馬氏鏈的蒙特卡羅方法，即Metropolis算法，並在最早的計算機上編程實現。Metropolis 算法是首個普適的采樣方法，並啟發了一系列 MCMC方法，所以人們把它視為隨機模擬技術騰飛的起點。Metropolis的這篇論文被收錄在《統計學中的重大突破》中，Metropolis算法也被遴選為二十世紀的十個最重要的算法之一。

我們接下來介紹的MCMC 算法是Metropolis算法的一個改進變種，即常用的 Metropolis-Hastings 算法。由上一節的例子和定理我們看到了，馬氏鏈的收斂性質主要由轉移矩陣P 決定, 所以基於馬氏鏈做采樣的關鍵問題是如何構造轉移矩陣P,使得平穩分布恰好是我們要的分布p(x)。如何能做到這一點呢？我們主要使用如下的定理。

如果非周期馬氏鏈的轉移矩陣\(T(z^*,z)\)和分布\(\pi(x)\)滿足$$\pi(i)T(i,j)=\pi(j)T(j,i)$$
則\(\pi(x)\)是馬氏鏈的平穩分布，上式被稱為細致平穩條件(detailed balance condition)。
滿足細致平穩條件就能收斂到平穩分布，平穩分布的定義

\begin{split}
&\pi(z)=\sum_{z^} T(z^,z)\
&\pi=\pi P
\end{split}

\begin{split}
&\pi(z)=\sum_{z^} T(z^,z)\pi(z^*)\
&\pi=\pi P
\end{split}

我們將細致平穩條件公式帶入到上式右邊，即可推出平穩分布定義

\[\sum_{z^*} T(z^*,z)\pi(z^*)=\sum_{z^*} T(z,z^*)\pi(z)=\pi(z)\sum_{z^*} p(z^*|z)=\pi(z) \]

1.Metropolis-Hastings

假設我們已經有一個轉移矩陣q(i,j)(表示從狀態i轉移到狀態j的概率，也可以寫為q(j|i)), 顯然，通常情況下$$p(i)q(i,j)\neq p(j)q(j,i)$$
也就是細致平穩條件不成立，所以p(x)不太可能是這個馬氏鏈的平穩分布。我們可否對馬氏鏈做一個改造，使得細致平穩條件成立呢？譬如，我們引入一個a(i,j) , 我們希望$$p(i)\underbrace{ q(i,j)a(i,j) }{q'(i,j)}=p(j)\underbrace{ q(j,i)a(j,i) }{q'(j,i)}$$
取什么樣的 以上等式能成立呢？最簡單的，按照對稱性，我們可以取$$a(i,j)=p(j)q(j,i)\qquad a(j,i)=p(i)q(i,j)$$
於是細致平穩條件就成立了，我們把原來具有轉移矩陣Q的一個很普通的馬氏鏈，改造為了具有轉移矩陣Q'的馬氏鏈，而Q'恰好滿足細致平穩條件，由此馬氏鏈的平穩分布就是p(x)!

在改造Q的過程中引入的a(i,j)稱為接受率，物理意義可以理解為在原來的馬氏鏈上，從狀態i以q(i,j)的概率轉跳轉到狀態j的時候，我們以a(i,j)的概率接受這個轉移，於是得到新的馬氏鏈Q'的轉移概率為q(i,j)a(i,j)。

假設我們已經有一個轉移矩陣Q(對應元素為q(i,j))，把以上的過程整理一下，我們就得到了如下的用於采樣概率分布的算法

• basic Metropolis
• 假設：直接從分布p(z)中采樣困難
• 問題：從分布p(z)中采樣
• 算法：

1. 首先初始化馬氏鏈的初始狀態\(X_0=x_0\)
2. 第t時刻狀態為\(x_t\)，對t+1時刻采樣\(y\sim p(x|x_t)\)
3. 從均勻分布采樣\(u\sim Uniform[0,1]\)
4. 如果\(u<a(x_t,y)=p(y)q(x_t|y)\)則接受轉移\(x_t\sim y\)，即\(X_{t+1}=y\)
5. 否則不接受轉移\(X_{t+1}=x_t\)
6. 跳轉第2步

• 缺陷：馬氏鏈在轉移的過程中的接受率a(i,j)可能偏小，這樣采樣過程中馬氏鏈容易原地踏步，拒絕大量的跳轉，這使得馬氏鏈遍歷所有的狀態空間要花費太長的時間，收斂到平穩分布p(x)的速度太慢。

有沒有辦法提升一些接受率呢?
假設a(i,j)=0.1 a(j,i)=0.2 , 此時滿足細致平穩條件，於是$$p(i)q(i,j)\times 0.1=p(j)q(j,i)\times 0.2$$
上式兩邊擴大5倍，我們改寫為$$p(i)q(i,j)\times 0.5=p(j)q(j,i)\times 1$$

看，我們提高了接受率，而細致平穩條件並沒有打破！這啟發我們可以把細致平穩條件式中的a(i,j)和a(j,i)同比例放大，使得兩數中最大的一個放大到1，這樣我們就提高了采樣中的跳轉接受率。所以我們可以取

\[a(i,j)=min\{\frac{p(j)q(j,i)}{p(i)q(i,j)},1\} \]

於是，經過對上述 basic Metropolis 采樣算法中接受率的微小改造，我們就得到了如下教科書中最常見的 Metropolis-Hastings 算法。

• Metropolis-Hastings
• 假設：直接從分布p(z)中采樣困難
• 問題：從分布p(z)中采樣
• 算法：

1. 首先初始化馬氏鏈的初始狀態\(X_0=x_0\)
2. 第t時刻狀態為\(x_t\)，對t+1時刻采樣\(y\sim p(x|x_t)\)
3. 從均勻分布采樣\(u\sim Uniform[0,1]\)
4. 如果\(u<a(x_t,y)=min\{\frac{p(y)q(x_t|y)}{p(x_t)q(y|x_t)},1\}\)則接受轉移\(x_t\sim y\)，即\(X_{t+1}=y\)
5. 否則不接受轉移\(X_{t+1}=x_t\)
6. 跳轉第2步
   對於分布p(x) ,我們構造轉移矩陣Q使其滿足細致平穩條件。

2.Gibbs Sampling

對於高維的情形，由於接受率a的存在(通常a<1), 以上 Metropolis-Hastings 算法的效率不夠高。能否找到一個轉移矩陣Q使得接受率a=1呢？我們先看看二維的情形，假設有一個概率分布p(x,y) ，考察x坐標相同的兩個點\(A(x_1,y_1)，A(x_1,y_2)\)我們發現

\begin{split}
& p(x_1,y_1)p(y_2|x_1)=p(x_1)p(y_1|x_1)p(y_2|x_1) \
& p(x_1,y_2)p(y_1|x_1)=p(x_1)p(y_2|x_1)p(y_1|x_1)
\end{split}

所以得到

\begin{split}
& p(x_1,y_1)p(y_2|x_1)=p(x_1,y_2)p(y_1|x_1) \qquad (1) \
& p(A)p(y_2|x_1)=p(B)p(y_1|x_1)
\end{split}

基於以上等式，我們發現，在\(x=x_1\)這條平行於y軸的直線上，如果使用條件分布\(p(y|x_1)\)做為任何兩個點之間的轉移概率，那么任何兩個點之間的轉移滿足細致平穩條件。同樣的，如果我們在\(y=y_1\)這條直線上任意取兩個點 ，也有如下等式

\[p(A)p(x2|y1)=p(C)p(x1|y1) \]

於是我們可以如下構造平面上任意兩點之間的轉移概率矩陣Q

\begin{split}
& Q(A,B)=p(y_B|x_1)\qquad 如果x_A=x_B=x_1\
& Q(A,C)=p(x_C|y_1)\qquad 如果y_A=y_C=y_1\
& Q(A,D)=0\qquad
\end{split}

有了如上的轉移矩陣Q, 我們很容易驗證對平面上任意兩點 , 滿足細致平穩條件

\[p(X)Q(X,Y)=p(Y)Q(Y,X) \]

於是這個二維空間上的馬氏鏈將收斂到平穩分布 。以上的過程我們很容易推廣到高維的情形，對於(1)式，如果\(x_1\)變為多維情形\(\vec x_1\)，可以看出推導過程不變，所以細致平穩條件同樣是成立的

\[p(\vec x_1,y_1)p(y_2|\vec x_1)=p(\vec x_1,y_2)p(y_1|\vec x_1) \]

此時轉移矩陣Q由條件分布\(p(y|\vec x_1)\)定義。上式只是說明了一根坐標軸的情形，和二維情形類似，很容易驗證對所有坐標軸都有類似的結論。所以n維空間中對於概率分布\(p(x_1,x_2,...,x_n)\)可以如下定義轉移矩陣

\begin{split}
& Q(A,B)=p(x_i|x_1,...,x_{i-1},x_{i+1}...,x_n)\qquad 如果沿着x_i這根坐標軸做轉移的時候\
& Q(A,D)=0\qquad others
\end{split}

• Gibbs Sampling
• 假設：直接從分布p(z)中采樣困難
• 問題：從分布p(z)中采樣
• 算法：

1. 首先初始化馬氏鏈的初始狀態\(X_0=x_0\)
2. 第t時刻狀態為\(x_t\)，對t+1時刻采樣
3. \(x_1^{(t+1)}\sim p(x_1|x_2^{(t)},x_3^{(t)},...,x_n^{(t)})\)
4. \(x_2^{(t+1)}\sim p(x_1|x_1^{(t+1)},x_3^{(t)},...,x_n^{(t)})\)
5. ...
6. \(x_j^{(t+1)}\sim p(x_1|x_1^{(t+1)},x_{j-1}^{(t+1)},...,x_{j+1}^{(t)},x_n^{(t)})\)
7. ...
8. \(x_n^{(t+1)}\sim p(x_1|x_1^{(t+1)},x_2^{(t)},...,x_{n-1}^{(t+1)})\)
9. 跳轉第3步

以上算法收斂后，得到的就是概率分布\(p(x_1,x_2,..,x_n)\)的樣本，當然這些樣本並不獨立，但是我們此處要求的是采樣得到的樣本符合給定的概率分布，並不要求獨立。同樣的，在以上算法中，坐標軸輪換采樣不是必須的，可以在坐標軸輪換中引入隨機性，這時候轉移矩陣 中任何兩個點的轉移概率中就會包含坐標軸選擇的概率，而在通常的Gibbs Sampling算法中，坐標軸輪換是一個確定性的過程，也就是在給定時刻，在一根固定的坐標軸上轉移的概率是1。

應用案例

1.LDA

參考資料

1. http://www.52nlp.cn/lda-math-mcmc-和-gibbs-sampling2
2. PRML, chapter 11
3. http://blog.quantitations.com/inference/2012/11/24/rejection-sampling-proof/
4. http://thexbar.me/2014/11/07/reject-sample/
5. Notes on Pattern Recognition and Machine Learning (Jian Xiao)
6. Pattern Recognition And Machine Learning 讀書會, chapter 11

Written with StackEdit.