MCMC 方法是從特定分布中生成隨機(jī)數(shù)的方法.? 不同于之前說過的逆變換采樣,? MCMC 方法可以適用于多維度分布,? 并且不存在過多額外的數(shù)學(xué)過程.

Markov Chain Monte Carlo 方法,? 顧名思義是一種依賴?Markov Chain (馬爾科夫鏈) 的隨機(jī)采樣方法,? 所以首先有必要來認(rèn)識一下馬爾科夫鏈.

馬爾科夫鏈?zhǔn)侵笍囊粋€(gè)狀態(tài)轉(zhuǎn)移到另一個(gè)狀態(tài)的隨機(jī)過程,? 并且轉(zhuǎn)移概率只與當(dāng)前狀態(tài)有關(guān),? 不與之前狀態(tài)有關(guān).? 對于這篇專欄來說了解到這樣子就可以了,? 詳細(xì)可以參考wiki.

至于 Monte Cario (蒙特卡洛) 涉及到積分的計(jì)算,? 然而這里不會計(jì)算積分所以就忽略掉了.

關(guān)于符號,? 這篇專欄里統(tǒng)一使用偏向數(shù)學(xué)風(fēng)格的符合記述,? 所以向量默認(rèn)列向量;??P(x=x?) 表示隨即狀態(tài) x?為 x? 的概率,? 簡記為 P(x?);? P(x=x?|y=y?) 表示?x?為 y? 時(shí) x 為 x?的概率,? 簡記為 P(x?|y?);? P(x=x?,y=y?) 表示 x 為 x? 并且 y 為 y? 的概率,? 簡記為 P(x?,y?);?etc.

隨機(jī)狀態(tài)的轉(zhuǎn)移和分布

因?yàn)轳R爾科夫鏈里狀態(tài)轉(zhuǎn)移的概率只與當(dāng)前狀態(tài)有關(guān),? 那么狀態(tài) i?轉(zhuǎn)移到狀態(tài) j 的概率可以表示為?,? 可以使用 轉(zhuǎn)移矩陣(transition matrix) 表示?.

使用 D??表示在隨機(jī)狀態(tài)在時(shí)刻 t 的分布(distibution),? D?(i) 表示隨機(jī)狀態(tài)在時(shí)刻 t 取得狀態(tài) i 的概率,? 那么下一時(shí)刻的分布 D??? 可以由 ?得出,? 或者寫為向量形式:?.

當(dāng)有一個(gè)特定的分布?π,? 經(jīng)過狀態(tài)轉(zhuǎn)移后仍然符合這個(gè)分布,? 即?,? 那么則這個(gè)分布稱為平穩(wěn)分布(stationary distribution).? 以向量形式來說為?,? 即?π 為轉(zhuǎn)移矩陣 T 的特征向量,? 并且特征值為1.

如果轉(zhuǎn)移矩陣 T 不可約并且是非周期的,? 那么以任意分布 D? 開始,? 執(zhí)行無限次轉(zhuǎn)移后,? 都會得到平穩(wěn)分布,? 即?,? 這個(gè)結(jié)論由 Perron–Frobenius 理論給出.??下面寫一個(gè)小程序?qū)@個(gè)結(jié)論進(jìn)行粗略的驗(yàn)證:

可以看到只要 T 不改變,? 無論是什么起始分布,? 最終都會收斂到 [0.625??0.3125?0.0625]?.

既然無論任何分布,? 施加多次轉(zhuǎn)移矩陣 T 都會使得分布趨向 T 對應(yīng)的平穩(wěn)分布?π,? 那么可以使用任意起始狀態(tài)的馬爾科夫鏈作為對平穩(wěn)分布的采樣.? 但需要注意的是,? 馬爾科夫鏈起始幾項(xiàng)不太符合平穩(wěn)分布,? 所以采樣需要在分布趨向平穩(wěn)后開始.? 下面的以上面的轉(zhuǎn)移矩陣作的例子

可以看到采樣分布比較符合平穩(wěn)分布.

MH 采樣

在一般情況下,? 我們是需要對某特定分布進(jìn)行采樣,? 而不是從轉(zhuǎn)移里得到分布,? 所以有必要從特定分布里求出轉(zhuǎn)移.

當(dāng)分布 π 為轉(zhuǎn)移?T 的平穩(wěn)分布時(shí),? 某狀態(tài) i 的概率密度流出量等于流入量,? 即滿足等式?,? 滿足這個(gè)等式的一個(gè)充分條件為 ,? 后一條件被稱為細(xì)致平衡(detailed balance)條件.

對于特定分布?π,? 選取任意轉(zhuǎn)移?Q,? 則存在?α 使得? 成立.? 不難知道,???是一種方便的取法.? 于是轉(zhuǎn)移??的平穩(wěn)分布即為?π.

使用接受-拒絕(acceptance-rejection)采樣的思路,? 生成相應(yīng)的馬爾科夫鏈步驟如下:? 1) 選取隨機(jī)起始狀態(tài) x?;? 2) 進(jìn)入循環(huán) t = 1, 2, ... :? 2.a) 從條件分布 Q(x'|x???) 中得到候選狀態(tài)?x';? 2.b) 計(jì)算轉(zhuǎn)移接受概率?α = Q(x???|x')?π(x');? 2.c) 選取分布在區(qū)間 [0,1) 的均勻隨機(jī)數(shù) s;? 2.d) 如果 s < α, 那么接受轉(zhuǎn)移?x? = x',? 否則不接受轉(zhuǎn)移?x? = x???.? 3) 最后截?cái)嗥鹗歼€沒符合分布 π 的鏈,? 得到可以代表分布?π 的采樣集.

[雖然這里可以插入一個(gè)例子,? 但是我真的很懶]

上述算法有兩個(gè)很嚴(yán)重的問題:? 看到接受概率 α,? 在離散情況時(shí)這是兩個(gè)小于 1 的數(shù)字相乘,? 這意味著接受率很小,? 從而導(dǎo)致需要很長的鏈才能得到足夠符合分布的采樣.? 另外,? 當(dāng)狀態(tài)不是離散而是連續(xù)時(shí),? 描述狀態(tài)分布由向量變?yōu)楹瘮?shù),? 即概率密度函數(shù)(probability density function),? 根據(jù) pdf 的定義,? pdf 的值允許大于1,? 所以接受概率 α 可能大于1,? 大于 1 的概率是沒有意義的.

為此,? 引入變量 b?使得?,??可以看到無論 b?取何值,? 細(xì)致平衡條件都是成立的.? 為了解決在離散時(shí) α 過小以及在連續(xù)時(shí) α 可能無意義兩個(gè)問題,? 選取 b?使得 ,? 不難得出 ,? 又記 ,? 那么 β 即是修改后的轉(zhuǎn)移接受概率.? 這種算法被稱為 Metropolis-Hastings 算法.

其他 MCMC 方法都可以看作是 MH 算法的特例.? 如果轉(zhuǎn)移 Q 的對稱的,? 即?,? MH 算法退化為 Metropolis 算法.

關(guān)于 Q 的選取一般有兩種:? 1) 直接選取一個(gè)與目標(biāo)分布接近且容易生成的分布;? 2) 以當(dāng)前狀態(tài)為中心向附近跳躍一段距離 (正態(tài)分布是好文明).

下面以第二種選取 Q 的方法作為例子:

可以看到生成的馬爾科夫鏈分布是非常接近目標(biāo)分布的,? 并且轉(zhuǎn)移接受率有 0.64.

不難知道,? MH 采樣極度依賴于轉(zhuǎn)移 Q,? 如果 Q 選取不當(dāng),? 接受轉(zhuǎn)移率降低會導(dǎo)致:? 1) 馬爾科夫鏈需要拋棄起始極長一段;? 2) 生成的隨機(jī)采樣質(zhì)量差;? 3) 計(jì)算成本增大.

選取合適的轉(zhuǎn)移最簡單就是觀察目標(biāo)分布的特征,? 然后選取一個(gè)接近且容易計(jì)算的分布.? 以上面展示的分布為例,? 轉(zhuǎn)移取 "均值=0, 標(biāo)準(zhǔn)差=7.5"的正態(tài)分布的右半部分,? 轉(zhuǎn)移接受率可以達(dá)到 0.47.? 對于 MH 算法來說轉(zhuǎn)移接受率大于 0.2 ~ 0.45 就已經(jīng)可以使用了,? 具體還是要看應(yīng)用場景.

另外還有一個(gè)問題,? 跟離散情況一樣,? 馬爾科夫鏈開頭幾項(xiàng)認(rèn)為是不符合目標(biāo)分布,? 所以需要拋棄的,? 但是究竟需要拋棄多長的鏈?? 可以同時(shí)運(yùn)行許多許多不同起始狀態(tài)的馬爾科夫鏈,? 抽取同一時(shí)刻的項(xiàng)跟目標(biāo)分布進(jìn)行對比,? 這樣可以大概地確定需要拋棄的數(shù)目.

高維采樣

需要注意的是,? "狀態(tài)" 并不是一個(gè)限制在一維的數(shù)字,? 實(shí)際上就像數(shù)學(xué)上定義 "集合" 為一堆東西組成的對象,? 并沒有對 "東西" 作定義,? 而 "狀態(tài)" 可以理解為集合里的元素,? 所以狀態(tài)實(shí)際上也是任意的.? 以下篇幅主要討論一下二維連續(xù)分布,? 更高維度的也可以作類比.

下面用 MH 采樣對分布? 進(jìn)行采樣,? 其中 N 是歸一化因子,? 眾所周知?dú)w一化因子不重要好吧.

可以看到采樣的分布比較符合目標(biāo)分布.? 但是有一個(gè)問題,? 轉(zhuǎn)移接受率只有 0.27,? 實(shí)際上這個(gè)數(shù)字對于高維的?MH 采樣來說已經(jīng)算不錯的了,? 普遍來說轉(zhuǎn)移接受率都在 0.1 左右.? 對此,? 這里提出 Gibbs 采樣.

在二維空間里進(jìn)行采樣,? 采樣步驟可以分解為:? 1) 先對 x 分量進(jìn)行采樣,? 2) 在取得 x 后再對 y 進(jìn)行采樣,? 即是?.? 但?x, y 的采樣順序是可以對換的,? 即有 .? 實(shí)際上聯(lián)合上面兩式可以得到貝葉斯定理:? .

考慮只有一個(gè)坐標(biāo)不相同的兩點(diǎn),? 比如說 x 不一樣:? (x?,y) 和 (x?,y),? 根據(jù)上面提出的采樣分解步驟,? 不難知道等式??成立,? 對比細(xì)致平衡條件,? 可以知道在 x 方向上轉(zhuǎn)移矩陣可以為 ,? 類似地 y 方向也有相同的結(jié)論:??.

于是 Gibbs 采樣步驟為:? (起始狀態(tài)和循環(huán)什么的就不寫了)? 1) 已知狀態(tài)?(x?,y?);? 2) 在 x 方向上根據(jù) Pr(x|y?) 采樣得到 x???;? 3) 在 y 方向上根據(jù) Pr(y|x???) 采樣得到 y???;? 4) 得到下一個(gè)狀態(tài)?(x???,y???).

在有些地方見到的 Gibbs 采樣在細(xì)節(jié)上會有不一樣的地方,? 比如說有些地方會把中間產(chǎn)生的狀態(tài) (x???,y?) 也當(dāng)作一個(gè)合理的采樣;? 還有些地方不是坐標(biāo)軸輪換采樣,? 而是隨機(jī)選取一個(gè)坐標(biāo)軸進(jìn)行采樣.? 這些變化都是沒問題的,? 完全ok.

以上面的分布??為例,? 可以求得條件分布: ,? 其中 ?分別是與 x 無關(guān)的歸一化因子,? 可以看到 x 對于 y 的條件分布是 "均值為 ,?標(biāo)準(zhǔn)差為 " 的正態(tài)分布,? y 對于 x 的分布也有相同的結(jié)論.? 于是得到 Gibbs 采樣的程序:

可以看到在 Gibbs 采樣里,? 因?yàn)榻邮?拒絕環(huán)節(jié)被取消了,? 生成的樣本比 MH 采樣更加均勻.? 實(shí)際上如果把中途生成的狀態(tài)?(x???,y?) 也當(dāng)作采樣,? 就可以看到狀態(tài)是如何轉(zhuǎn)移的:

Gbbis 采樣相比與 MH 采樣,? 優(yōu)勢就是以 100% 概率接受狀態(tài)轉(zhuǎn)移,? 并且不需要知道多個(gè)變量之間的聯(lián)合分布,? 只需要知道單個(gè)變量對于其他變量的條件分布,? 在一些情況下聯(lián)合分布比條件分布難求得多.

不過 Gibbs 采樣也不是萬能的,??如果條件分布也是復(fù)雜分布的話,? 雖然可以使用簡單的接受-拒絕采樣,? 但是從成本上來說說不定比 MH 采樣更差.? 對于條件分布極難甚至不能求出的情況,? Gbbis 采樣就完全不適用,? 這時(shí)候可以采用更高級的采樣方法,? 這里舉幾個(gè)例子:? 自適應(yīng) MH 采樣,??Hamiltonian 方法,??Langevin 規(guī)則,? 切片采樣,? 可逆跳躍采樣.

采樣方法多種多樣,? 不過還是那句話

所以更高級的采樣方法我也懶得去翻了,? 這篇專欄只是學(xué)量子力學(xué)途中遇上一點(diǎn)難題,? 才順便拉出來的.

日常推銷瑟圖群:??274767696

封面pid:?66081482