從隨機過程到馬爾科夫鏈蒙特卡洛方法

時間 2019-11-10

標籤隨機過程蒙特卡洛方法简体版

原文原文鏈接

1. Introduction

第一次接觸到 Markov Chain Monte Carlo (MCMC) 是在 theano 的 deep learning tutorial 裏面講解到的 RBM 用到了 Gibbs sampling，當時由於要趕着作項目，雖然一頭霧水，可是也沒沒有時間仔細看。趁目前比較悠閒，把 machine learning 裏面的 sampling methods 理一理，發現內容還真很多，有些知識本人也是隻知其一;不知其二，因此這篇博客不可能面面俱到詳細講解全部的 sampling methods，而是着重講一下這個號稱二十世紀 top 10 之一的算法—— Markov chain Monte Carlo。在介紹 MCMC 以前，咱們首先了解一下 MCMC 的 Motivation 和在它以前用到的方法。本人也是初學者，錯誤在所不免，歡迎一塊兒交流。python

這篇博客從零開始，應該均可以看懂，主要內容包括：git

隨機採樣
拒絕採樣
重要性採樣
Metropolis-Hastings Algorithm
Gibbs Sampling

2. Sampling

咱們知道，計算機自己是沒法產生真正的隨機數的，可是能夠根據必定的算法產生僞隨機數（pseudo-random numbers)。最古老最簡單的莫過於 Linear congruential generator：算法

式子中的 a 和 c 是一些數學知識推導出的合適的常數。可是咱們看到，這種算法產生的下一個隨機數徹底依賴如今的隨機數的大小，並且當你的隨機數序列足夠大的時候，隨機數將出現重複子序列的狀況。固然，理論發展到今天，有不少更加先進的隨機數產生算法出現，好比 python 數值運算庫 numpy 用的是 Mersenne Twister 等。可是無論算法如何發展，這些都不是本質上的隨機數，用馮諾依曼的一句話說就是：網頁爬蟲

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

要檢查一個序列是不是真正的隨機序列，能夠計算這個序列的 entropy 或者用壓縮算法計算該序列的冗餘。app

OK，根據上面的算法如今咱們有了均勻分佈的隨機數，可是如何產生知足其餘分佈（好比高斯分佈）下的隨機數呢？一種可選的簡單的方法是 Inverse transform sampling，有時候也叫Smirnov transform。拿高斯分佈舉例子，它的原理是利用高斯分佈的累積分佈函數（CDF，cumulative distribution function）來處理，過程以下圖：dom

假如在 y 軸上產生（0,1）之間的均勻分佈的隨機數，水平向右投影到高斯累計分佈函數上，而後垂直向下投影到 x 軸，獲得的就是高斯分佈。可見高斯分佈的隨機數實際就是均勻分佈隨機數在高斯分佈的 CDF 函數下的逆映射。固然，在實際操做中，更有效的計算方法有 Box–Muller_transform (an efficient polar form)，Ziggurat algorithm 等，這些方法 tricky and faster，沒有深刻了解，這裏也很少說了。ide

3. Motivation

MCMC 可解決高維空間裏的積分和優化問題：函數

上面一個例子簡單講了利用高斯分佈的 CDF 能夠產生高斯隨機數，可是有時候咱們遇到一些分佈的 CDF 計算不出來（沒法用公式表示），隨機數如何產生？優化
遇到某些沒法直接求積分的函數，如 e^{x^2}，在計算機裏面如何求積分？
如何對一個分佈進行高效快速的模擬，以便於抽樣？
如何在可行域很大(or large number of possible configurations)時有效找到最優解——RBM 優化目標函數中的問題。
如何在衆多模型中快速找到更好的模型——MDL, BIC, AIC 模型選擇問題。

3.1 The Monte Carlo principle

實際上，Monte Carlo 抽樣基於這樣的思想：假設玩一局牌的贏的機率只取決於你抽到的牌，若是用窮舉的方法則有 52! 種狀況，計算複雜度太大。而現實中的作法是先玩幾局試試，統計贏的機率，若是你不太確信這個機率，你能夠儘量多玩幾局，當你玩的次數很大的時候，獲得的機率就很是接近真實機率了。

上述方法能夠估算隨機事件的機率，而用 Monte Carlo 抽樣計算隨即變量的指望值是接下來內容的重點：X 表示隨即變量，服從機率分佈 p(x), 那麼要計算 f(x) 的指望，只須要咱們不停從 p(x) 中抽樣

當抽樣次數足夠的時候，就很是接近真實值了：

Monte Carlo 抽樣的方法還有一個重要的好處是：估計值的精度與 x 的維度無關（雖然維度越高，可是每次抽樣得到的信息也越多），而是與抽樣次數有關。在實際問題裏面抽樣二十次左右就能達到比較好的精度。

可是，當咱們實際動手的時候，就會發現一個問題——如何從分佈 p(x) 中抽取隨機樣本。以前咱們說過，計算能夠產生均勻分佈的僞隨機數。顯然，第二小節產生高斯隨機數的抽樣方法只對少數特定的問題管用，對於通常狀況呢？

3.2 Rejection Sampling

Reject Sampling 實際採用的是一種迂迴( proposal distribution q(x) )的策略。既然 p(x) 太複雜在程序中無法直接採樣，那麼我設定一個程序可抽樣的分佈 q(x) 好比高斯分佈，而後按照必定的方法拒絕某些樣本，達到接近 p(x) 分佈的目的：

具體操做以下，設定一個方便抽樣的函數 q(x)，以及一個常量 k，使得 p(x) 總在 kq(x) 的下方。（參考上圖）

x 軸方向：從 q(x) 分佈抽樣獲得 a。(若是是高斯，就用以前說過的 tricky and faster 的算法更快）
y 軸方向：從均勻分佈（0, kq(a)) 中抽樣獲得 u。
若是恰好落到灰色區域： u > p(a), 拒絕，不然接受此次抽樣
重複以上過程

用均勻分佈拒絕抽樣來近似兩個高斯混合分佈的代碼以下：

rejectionsampling.py

# -*- coding=utf8 -*-

# Code from Chapter 14 of Machine Learning: An Algorithmic Perspective
# The basic rejection sampling algorithm

from pylab import *
from numpy import *

def qsample():
    return random.rand()*4.

def p(x):
    return 0.3*exp(-(x-0.3)**2) + 0.7* exp(-(x-2.)**2/0.3) 

def rejection(nsamples):
    
    M = 0.72#0.8
    samples = zeros(nsamples,dtype=float)
    count = 0
    for i in range(nsamples):
        accept = False
        while not accept:
            x = qsample()
            u = random.rand()*M
            if u<p(x):
                accept = True
                samples[i] = x
            else: 
                count += 1
    print count   
    return samples

x = arange(0,4,0.01)
x2 = arange(-0.5,4.5,0.1)
realdata = 0.3*exp(-(x-0.3)**2) + 0.7* exp(-(x-2.)**2/0.3) 
box = ones(len(x2))*0.75#0.8
box[:5] = 0
box[-5:] = 0
plot(x,realdata,'k',lw=6)
plot(x2,box,'k--',lw=6)

import time
t0=time.time()
samples = rejection(10000)
t1=time.time()
print "Time ",t1-t0

hist(samples,15,normed=1,fc='k')
xlabel('x',fontsize=24)
ylabel('p(x)',fontsize=24)
axis([-0.5,4.5,0,1])
show()

View Code

在高維的狀況下，Rejection Sampling 會出現兩個問題，第一是合適的 q 分佈比較難以找到，第二是很難肯定一個合理的 k 值。這兩個問題會致使拒絕率很高，無用計算增長。

3.3 Importance Sampling

Importance Sampling 也是藉助了容易抽樣的分佈 q (proposal distribution)來解決這個問題，直接從公式出發：

其中，p(z) / q(z) 能夠看作 importance weight。咱們來考察一下上面的式子，p 和 f 是肯定的，咱們要肯定的是 q。要肯定一個什麼樣的分佈纔會讓採樣的效果比較好呢？直觀的感受是，樣本的方差越小指望收斂速率越快。好比一次採樣是 0, 一次採樣是 1000, 平均值是 500,這樣採樣效果不好，若是一次採樣是 499, 一次採樣是 501, 你說指望是 500,可信度還比較高。在上式中，咱們目標是 p×f/q 方差越小越好，因此 |p×f| 大的地方，proposal distribution q(z) 也應該大。舉個稍微極端的例子：

第一個圖表示 p 分佈，第二個圖的陰影區域 f = 1，非陰影區域 f = 0, 那麼一個良好的 q 分佈應該在左邊箭頭所指的區域有很高的分佈機率，由於在其餘區域的採樣計算實際上都是無效的。這代表 Importance Sampling 有可能比用原來的 p 分佈抽樣更加有效。

可是惋惜的是，在高維空間裏找到一個這樣合適的 q 很是難。即便有 Adaptive importance sampling 和 Sampling-Importance-Resampling(SIR) 的出現，要找到一個同時知足 easy to sample 而且 good approximations 的 proposal distribution, it is often impossible！

4. MCMC Algorithm

上面說了這麼多采樣方法，其實最終要突出的就是 MCMC 的過人之處。MCMC 的絕妙之處在於：經過穩態的 Markov Chain 進行轉移計算，等效於從 P(x) 分佈採樣。可是在瞭解 MCMC 具體算法以前，咱們還要熟悉 Markov Chain 是怎麼一回事。

4.1 Markov Chain

Markov Chain 體現的是狀態空間的轉換關係，下一個狀態只決定與當前的狀態(能夠聯想網頁爬蟲原理，根據當前頁面的超連接訪問下一個網頁)。以下圖：

這個狀態圖的轉換關係能夠用一個轉換矩陣 T 來表示：

舉一個例子，若是當前狀態爲 u(x) = (0.5, 0.2, 0.3), 那麼下一個矩陣的狀態就是 u(x)T = (0.18, 0.64, 0.18), 依照這個轉換矩陣一直轉換下去，最後的系統就趨近於一個穩定狀態 (0.22, 0.41, 0.37) (此處只保留了兩位有效數字)。而事實證實不管你從那個點出發，通過很長的 Markov Chain 以後都會聚集到這一點。

考慮通常的狀況，知足什麼條件下通過很長的 Markov Chain 迭代後系統分佈會趨近一個穩定分佈，即最後的 u(x) 等效於從目標分佈 p(x) 採樣。大概的條件以下（本身隨便總結的，可能有遺漏和錯誤）：

Irreducibility. 即圖是聯通的，T 矩陣不能被切豆腐同樣劃分紅小方塊，舉個例子，好比爬蟲爬不到內部局域網的網頁
Aperiodicity. 即圖中遍歷不會陷入到一個死圈裏，有些網站爲了防機器人，會專門設置這種陷阱
Detailed Balance，這是保證系統有穩態的一個重要條件，詳細說明見下面。

假設 p(x) 是最後的穩態，那麼 detailed balance 能夠用公式表示爲：

什麼意思呢？假設上面狀態圖 x1 有 0.22 元， x2 有 0.41 元，x3 有 0.37 元，那麼 0.22×1 表示 x1 須要給 x2 錢，以此類推，手動計算，能夠發現下一個狀態每一個人手中的錢都沒有變。值得說明的是，這裏體現了一個很重要的特性，那就是從一個高几率狀態 xi 向一個低機率狀態 x(i-1) 轉移的機率等於從這個低機率狀態向高几率狀態轉移的機率（reversible，至於要不要轉移又是另一回事)。固然，在上面一個例子中，狀況比較特殊，等號兩邊其實都是同一個東西。馬氏鏈的收斂性質主要由轉移矩陣決定, 因此基於馬氏鏈作採樣的關鍵問題是如何構造轉移矩陣,使得平穩分佈剛好是咱們要的分佈p(x)。可是考慮一維的狀況，假設 p(x) 是一維高斯分佈，x 是根據 markov chain 獲得的一個樣本，依照上面的等式，那麼咱們能夠根據轉移矩陣 T左和 T右（這裏實際是 proposal distribution）來獲得 p(xi) 和 p(x(i-1)) 的比率，進而按照必定的機率對這兩個樣本進行選擇。經過大量這樣的處理，獲得樣本就符合原始的 p(x) 分佈了。這就是 MH 算法的基本原理。

4.2 Metropolis-Hastings Algorithm

舉個例子，咱們要用 MH 算法對標準高斯分佈進行採樣，轉移函數(對稱)是方差爲 0.05 的高斯矩陣，上述算法過程以下：

選取一個隨機點 x0，做爲一個採樣點
以 x0 爲中心，以轉移函數爲分佈採起隨機點 x1
以算法中的 A 機率接受 x1, 不然接受 x0
重複第二步第三步

注意到高斯分佈是一個徑向基函數，上面算法畫波浪線的部分相等。

matlab 代碼以下：

n = 250000;
x = zeros(n, 1);
x(1) = 0.5;
for i = 1: n-1
    x_c = normrnd(x(1), 0.05);
    if rand < min(1, normpdf(x_c)/normpdf(x(i)))
        x(i+1) = x_c;
    else
        x(i+1) = x(i);
    end
end

MH 算法中的 proposal distribution q(x) 也是須要當心肯定的，詳細知識能夠查閱這篇介紹論文 (An introduction to MCMC for machine learning, Andrieu, Christophe). 能夠看到，這個算法和模擬退火算法的思想是很是類似的，可是在模擬退火算法過程當中，隨着時間的增長，接受值大的區域的機率愈來愈高，直到找到最高點。

4.3 Gibbs Sampling

Gibbs Sampling 其實是 MH 算法的一個變種。具體思路以下：假設在必定溫度下必定量的分子在容器裏作無規則的熱運動，如何統計系統的能量呢？一樣，咱們用 Monte Carlo 的思想進行統計計算。咱們假設全部的分子靜止在某一個時刻，這是初識狀態。固定其餘的分子，根據分子間的做用力對其中一個分子進行移動，也就是說在該分子以必定的機率移動到領域的某一個地方，移動完了以後再靜止。而後基於移動後的狀態對下一個分子進行一樣的移動操做...直到全部的分子移動完畢，那麼如今的狀態就是 Monte Carlo 採樣的第二個樣本。依照這樣的順序採樣下去，咱們對於這個系統就能計算一個統計意義上的能量了。從條件分佈的角度來看，算法過程以下：

整體來說，Gibbs Sampling 就是從條件機率中選擇一個變量（分子），而後對該變量（分子）進行採樣。當全部變量採樣完畢以後，就獲得了後面的一個狀態，從而完成了對系統配置的採樣。在 deep learning 的 RBM 中，gibbs 採樣是已知權重參數和一個 v 變量，經過採樣獲得 h。經過 h 採樣又能夠獲得另外一個 v ，如此交替採樣，就能夠逐漸收斂於聯合分佈了。

Gibbs.py (對高斯分佈進行 Gibbs 採樣)

# -*- coding=utf8 -*-

# Code from Chapter 14 of Machine Learning: An Algorithmic Perspective
# A simple Gibbs sampler

from pylab import *
from numpy import *

def pXgivenY(y,m1,m2,s1,s2):
    return random.normal(m1 + (y-m2)/s2,s1)

def pYgivenX(x,m1,m2,s1,s2):
    return random.normal(m2 + (x-m1)/s1,s2)

def gibbs(N=5000):
    k=20
    x0 = zeros(N,dtype=float)
    m1 = 10
    m2 = 20
    s1 = 2
    s2 = 3
    for i in range(N):
        y = random.rand(1)
        # 每次採樣須要迭代 k 次
        for j in range(k):
            x = pXgivenY(y,m1,m2,s1,s2)
            y = pYgivenX(x,m1,m2,s1,s2)
        x0[i] = x
    
    return x0

def f(x):
    return exp(-(x-10)**2/10)

# 畫圖
N=10000
s=gibbs(N)
x1 = arange(0,17,1)
hist(s,bins=x1,fc='k')
x1 = arange(0,17, 0.1)
px1 = zeros(len(x1))
for i in range(len(x1)):
    px1[i] = f(x1[i])
plot(x1, px1*N*10/sum(px1), color='k',linewidth=3)

show()