簡(jiǎn)單易學(xué)的機(jī)器學(xué)習(xí)算法—馬爾可夫鏈蒙特卡羅方法MCMC

對(duì)于一般的分布的采樣，在很多的編程語(yǔ)言中都有實(shí)現(xiàn)，如最基本的滿足均勻分布的隨機(jī)數(shù)，但是對(duì)于復(fù)雜的分布，要想對(duì)其采樣，卻沒(méi)有實(shí)現(xiàn)好的函數(shù)，在這里，可以使用馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo, MCMC)方法，其中Metropolis-Hastings采樣和Gibbs采樣是MCMC中使用較為廣泛的兩種形式。
MCMC的基礎(chǔ)理論為馬爾可夫過(guò)程，在MCMC算法中，為了在一個(gè)指定的分布上采樣，根據(jù)馬爾可夫過(guò)程，首先從任一狀態(tài)出發(fā)，模擬馬爾可夫過(guò)程，不斷進(jìn)行狀態(tài)轉(zhuǎn)移，最終收斂到平穩(wěn)分布。
一、馬爾可夫鏈
1、馬爾可夫鏈
設(shè)Xt表示隨機(jī)變量X在離散時(shí)間t時(shí)刻的取值。若該變量隨時(shí)間變化的轉(zhuǎn)移概率僅僅依賴于它的當(dāng)前取值，即

也就是說(shuō)狀態(tài)轉(zhuǎn)移的概率只依賴于前一個(gè)狀態(tài)。稱這個(gè)變量為馬爾可夫變量，其中，s0,s1,?,si,sj∈Ω為隨機(jī)變量X可能的狀態(tài)。這個(gè)性質(zhì)稱為馬爾可夫性質(zhì)，具有馬爾可夫性質(zhì)的隨機(jī)過(guò)程稱為馬爾可夫過(guò)程。

馬爾可夫鏈指的是在一段時(shí)間內(nèi)隨機(jī)變量X的取值序列(X0,X1,?,Xm)，它們滿足如上的馬爾可夫性質(zhì)。

2、轉(zhuǎn)移概率

馬爾可夫鏈?zhǔn)峭ㄟ^(guò)對(duì)應(yīng)的轉(zhuǎn)移概率定義的，轉(zhuǎn)移概率指的是隨機(jī)變量從一個(gè)時(shí)刻到下一個(gè)時(shí)刻，從狀態(tài)si轉(zhuǎn)移到另一個(gè)狀態(tài)sj的概率，即：

記表示隨機(jī)變量X在時(shí)刻t的取值為sk的概率，則隨機(jī)變量X在時(shí)刻t+1的取值為si的概率為：

假設(shè)狀態(tài)的數(shù)目為n，則有：

3、馬爾可夫鏈的平穩(wěn)分布

對(duì)于馬爾可夫鏈，需要注意以下的兩點(diǎn)：

1、周期性：即經(jīng)過(guò)有限次的狀態(tài)轉(zhuǎn)移，又回到了自身；
2、不可約：即兩個(gè)狀態(tài)之間相互轉(zhuǎn)移；
如果一個(gè)馬爾可夫過(guò)程既沒(méi)有周期性，又不可約，則稱為各態(tài)遍歷的。

對(duì)于一個(gè)各態(tài)遍歷的馬爾可夫過(guò)程，無(wú)論初始值π(0)取何值，隨著轉(zhuǎn)移次數(shù)的增多，隨機(jī)變量的取值分布最終都會(huì)收斂到唯一的平穩(wěn)分布π?，即：

且這個(gè)平穩(wěn)分布π?滿足：

其中，為轉(zhuǎn)移概率矩陣。

二、馬爾可夫鏈蒙特卡羅方法

1、基本思想

對(duì)于一個(gè)給定的概率分布P(X)，若是要得到其樣本，通過(guò)上述的馬爾可夫鏈的概念，我們可以構(gòu)造一個(gè)轉(zhuǎn)移矩陣為P的馬爾可夫鏈，使得該馬爾可夫鏈的平穩(wěn)分布為P(X)，這樣，無(wú)論其初始狀態(tài)為何值，假設(shè)記為x0，那么隨著馬爾科夫過(guò)程的轉(zhuǎn)移，得到了一系列的狀態(tài)值，如：x0,x1,x2,?,xn,xn+1,?,，如果這個(gè)馬爾可夫過(guò)程在第n步時(shí)已經(jīng)收斂，那么分布P(X)的樣本即為xn,xn+1,?。

2、細(xì)致平穩(wěn)條件

對(duì)于一個(gè)各態(tài)遍歷的馬爾可夫過(guò)程，若其轉(zhuǎn)移矩陣為P，分布為π(x)，若滿足：

則π(x)是馬爾可夫鏈的平穩(wěn)分布，上式稱為細(xì)致平穩(wěn)條件。

3、Metropolis采樣算法

Metropolis采樣算法是最基本的基于MCMC的采樣算法。

3.1、Metropolis采樣算法的基本原理

假設(shè)需要從目標(biāo)概率密度函數(shù)p(θ)中進(jìn)行采樣，同時(shí)，θ滿足?∞<θ<∞。Metropolis采樣算法根據(jù)馬爾可夫鏈去生成一個(gè)序列：

其中，θ(t)表示的是馬爾可夫鏈在第t代時(shí)的狀態(tài)。
在Metropolis采樣算法的過(guò)程中，首先初始化狀態(tài)值θ(1)，然后利用一個(gè)已知的分布生成一個(gè)新的候選狀態(tài)θ(?)，隨后根據(jù)一定的概率選擇接受這個(gè)新值，或者拒絕這個(gè)新值，在Metropolis采樣算法中，概率為：

這樣的過(guò)程一直持續(xù)到采樣過(guò)程的收斂，當(dāng)收斂以后，樣本θ(t)即為目標(biāo)分布p(θ)中的樣本。

3.2、Metropolis采樣算法的流程

基于以上的分析，可以總結(jié)出如下的Metropolis采樣算法的流程：

初始化時(shí)間t=1
設(shè)置u的值，并初始化初始狀態(tài)θ(t)=u
重復(fù)一下的過(guò)程：
令t=t+1
從已知分布中生成一個(gè)候選狀態(tài)θ(?)
計(jì)算接受的概率：
從均勻分布Uniform(0,1)生成一個(gè)隨機(jī)值a
如果a?α，接受新生成的值：θ(t)=θ(?)；否則：θ(t)=θ(t?1)
直到t=T
3.3、Metropolis算法的解釋

要證明Metropolis采樣算法的正確性，最重要的是要證明構(gòu)造的馬爾可夫過(guò)程滿足如上的細(xì)致平穩(wěn)條件，即：

對(duì)于上面所述的過(guò)程，分布為p(θ)，從狀態(tài)i轉(zhuǎn)移到狀態(tài)j的轉(zhuǎn)移概率為：

其中，Qi,j為上述已知的分布。

對(duì)于選擇該已知的分布，在Metropolis采樣算法中，要求該已知的分布必須是對(duì)稱的，即Qi,j=Qj,i，即

常用的符合對(duì)稱的分布主要有：正態(tài)分布，柯西分布以及均勻分布等。
接下來(lái)，需要證明在Metropolis采樣算法中構(gòu)造的馬爾可夫鏈滿足細(xì)致平穩(wěn)條件。

因此，通過(guò)以上的方法構(gòu)造出來(lái)的馬爾可夫鏈?zhǔn)菨M足細(xì)致平穩(wěn)條件的。

3.4、實(shí)驗(yàn)

假設(shè)需要從柯西分布中采樣數(shù)據(jù)，我們利用Metropolis采樣算法來(lái)生成樣本，其中，柯西分布的概率密度函數(shù)為：

那么，根據(jù)上述的Metropolis采樣算法的流程，接受概率α的值為：

代碼如下：

'''
Date:20160629
@author: zhaozhiyong
'''
import random
from scipy.stats import norm
import matplotlib.pyplot as plt

def cauchy(theta):
    y = 1.0 / (1.0 + theta ** 2)
    return y

T = 5000
sigma = 1
thetamin = -30
thetamax = 30
theta = [0.0] * (T+1)
theta[0] = random.uniform(thetamin, thetamax)

t = 0
while t < T:
    t = t + 1
    theta_star = norm.rvs(loc=theta[t - 1], scale=sigma, size=1, random_state=None)
    #print theta_star
    alpha = min(1, (cauchy(theta_star[0]) / cauchy(theta[t - 1])))

    u = random.uniform(0, 1)
    if u <= alpha:
        theta[t] = theta_star[0]
    else:
        theta[t] = theta[t - 1]

ax1 = plt.subplot(211)
ax2 = plt.subplot(212)
plt.sca(ax1)
plt.ylim(thetamin, thetamax)
plt.plot(range(T+1), theta, 'g-')
plt.sca(ax2)
num_bins = 50
plt.hist(theta, num_bins, normed=1, facecolor='red', alpha=0.5)
plt.show()數(shù)據(jù)分析師培訓(xùn)

實(shí)驗(yàn)的結(jié)果：

對(duì)于Metropolis采樣算法，其要求選定的分布必須是對(duì)稱的。