【机器学习基础】数学推导+纯Python实现机器学习算法21：马尔可夫链蒙特卡洛...
Python機(jī)器學(xué)習(xí)算法實(shí)現(xiàn)
Author：louwill
Machine Learning Lab
? ? ?
? ? ?蒙特卡洛（Monte Carlo，MC）方法作為一種統(tǒng)計(jì)模擬和近似計(jì)算方法，是一種通過對(duì)概率模型隨機(jī)抽樣進(jìn)行近似數(shù)值計(jì)算的方法。馬爾可夫鏈（Markov Chain，MC）則是一種具備馬爾可夫性的隨機(jī)序列。將二者結(jié)合起來(lái)便有了馬爾可夫鏈蒙特卡洛方法（Markov Chain Monte Carlo，MCMC），即是以構(gòu)造馬爾可夫鏈為概率模型的蒙特卡洛方法。
? ? ?限于篇幅，本文不對(duì)MCMC的前置知識(shí)進(jìn)行詳細(xì)介紹，關(guān)于蒙特卡洛方法和馬爾可夫鏈的大量?jī)?nèi)容，請(qǐng)各位讀者自行查閱相關(guān)材料進(jìn)行學(xué)習(xí)。本文聚焦于MCMC方法本身原理和常用實(shí)現(xiàn)方法。
MCMC簡(jiǎn)介
? ? ?一般來(lái)說(shuō)，對(duì)目標(biāo)概率模型進(jìn)行隨機(jī)抽樣能夠幫助我們得到該分布的近似數(shù)值解。但如果隨機(jī)變量的多元的，或者所要抽樣的概率密度函數(shù)形式是復(fù)雜的非標(biāo)準(zhǔn)式時(shí)，直接應(yīng)用蒙特卡洛方法就會(huì)很困難。
? ? ?MCMC方法的基本思路是：在隨機(jī)變量的狀態(tài)空間上定義一個(gè)馬爾可夫鏈，使其平穩(wěn)分布就是我們要抽樣的目標(biāo)分布。然后基于該馬爾可夫鏈進(jìn)行隨機(jī)游走產(chǎn)生對(duì)應(yīng)的樣本序列進(jìn)行數(shù)值計(jì)算。當(dāng)時(shí)間足夠長(zhǎng)時(shí)，抽樣所得到的分布就會(huì)趨近于平穩(wěn)分布，基于該馬爾可夫鏈的抽樣結(jié)果就是目標(biāo)概率分布的抽樣結(jié)果。對(duì)抽樣結(jié)果的函數(shù)均值計(jì)算就是目標(biāo)數(shù)學(xué)期望值。
我們將上述流程梳理一下，完整的MCMC方法包括以下3個(gè)步驟：
在隨機(jī)變量的狀態(tài)空間上定義一個(gè)滿足遍歷定理的馬爾可夫鏈，使其平穩(wěn)分布為目標(biāo)分布；
從狀態(tài)空間種某一點(diǎn)出發(fā)，在構(gòu)造的馬爾可夫鏈上進(jìn)行隨機(jī)游走，可以得到樣本序列；
根據(jù)馬爾可夫鏈的遍歷定理，確定正整數(shù)和，可得樣本集合，最后可得函數(shù)的遍歷均值：
? ? ?所以MCMC的關(guān)鍵問題在于如何構(gòu)造滿足條件的馬爾可夫鏈。常用的MCMC構(gòu)建算法包括Metropolis-Hasting算法和Gibbs抽樣。
Metropolis-Hasting算法
? ? ?Metropolis-Hasting算法也可以簡(jiǎn)稱為M-H采樣，該算法由Metropolis提出后經(jīng)Hasting改進(jìn)，故而得名。假設(shè)目標(biāo)抽樣分布為，M-H算法采用的狀態(tài)轉(zhuǎn)移核為，則其馬爾可夫鏈可定義為：
? ? ?其中和分別為建議分布和接受分布。建議分布可以是另一個(gè)馬爾可夫鏈的轉(zhuǎn)移核，其形式也可能有多種。包括對(duì)稱形式和獨(dú)立抽樣形式。假設(shè)建議分布是對(duì)稱的，即對(duì)任意的和，有：
? ? ?接受分布形式如下：
? ? ?轉(zhuǎn)移核為的馬爾可夫鏈的隨機(jī)游走過程如下：如果在時(shí)刻處于狀態(tài)，即有，則先按照建議分布抽樣產(chǎn)生一個(gè)候選狀態(tài)，然后按照接受分布抽樣決定是否接受狀態(tài)。以概率接受，以概率拒絕。完整的Metropolis-Hasting算法流程如下：
在狀態(tài)空間上任意選擇一個(gè)初始值；
對(duì)遍歷執(zhí)行
設(shè)狀態(tài)，按照建議分布抽取一個(gè)候選狀態(tài)。
計(jì)算接受概率
從區(qū)間上按均勻分布隨機(jī)抽取一個(gè)數(shù)，若，則狀態(tài)；否則。
最后得到樣本集合，計(jì)算：
? ? ?借助Python的高級(jí)計(jì)算庫(kù)scipy，下面給出Metropolis-Hasting算法的基本實(shí)現(xiàn)過程。假設(shè)目標(biāo)平穩(wěn)分布是一個(gè)正態(tài)分布，基于M-H的采樣過程如下。
from scipy.stats import norm import random import matplotlib.pyplot as plt# 定義平穩(wěn)分布 def smooth_dist(theta):y = norm.pdf(theta, loc=3, scale=2)return yT = 10000 pi = [0 for i in range(T)] sigma = 1 # 設(shè)置初始值 t = 0 # 遍歷執(zhí)行 while t < T-1:t = t + 1# 狀態(tài)轉(zhuǎn)移進(jìn)行隨機(jī)抽樣pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None) # 計(jì)算接受概率alpha = min(1, (smooth_dist(pi_star[0]) / smooth_dist(pi[t - 1]))) # 從均勻分布中隨機(jī)抽取一個(gè)數(shù)uu = random.uniform(0, 1)# 拒絕-接受采樣if u < alpha:pi[t] = pi_star[0]else:pi[t] = pi[t - 1]# 繪制采樣分布 plt.scatter(pi, norm.pdf(pi, loc=3, scale=2), label='Target Distribution') num_bins = 100 plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.6, label='Samples Distribution') plt.legend() plt.show();經(jīng)過10000次遍歷迭代后，采樣效果如下圖所示。
Gibbs抽樣
? ? ?相較于M-H抽樣，Gibbs抽樣是目前更常用的MCMC抽樣算法，Gibbs抽樣可以視為特殊的M-H抽樣方法。Gibbs抽樣適用于多元隨機(jī)變量聯(lián)合分布的抽樣和估計(jì)，其基本思路是從聯(lián)合概率分布種定義滿條件概率分布，依次對(duì)滿條件概率分布進(jìn)行抽樣得到目標(biāo)樣本序列。
? ? ?這里簡(jiǎn)單提一下滿條件分布。假設(shè)MCMC的目標(biāo)分布為多元聯(lián)合概率分布，如果條件概率分布中所有個(gè)變量全部出現(xiàn)，其中，，這種條件概率分布即為滿條件分布。
? ? ?假設(shè)在第步得到樣本，在第步先對(duì)第一個(gè)隨機(jī)變量按照如下滿條件分布進(jìn)行隨機(jī)抽樣得到:
之后依次對(duì)第個(gè)變量按照滿條件概率分布進(jìn)行抽樣，最后可得整體樣本。
? ? ?Gibbs抽樣可視為單分量M-H抽樣的特殊情形。即Gibbs抽樣對(duì)每次抽樣結(jié)果都以1的概率進(jìn)行接受，從不拒絕，這跟M-H采樣有本質(zhì)區(qū)別。Gibbs抽樣的完整流程如下：
給定多隨機(jī)變量初始值。
對(duì)遍歷執(zhí)行：假設(shè)第步得到樣本，則第次迭代進(jìn)行如下步驟：
由滿條件分布抽取得到
由滿條件分布抽取
由滿條件分布抽取
最后得到樣本集合，計(jì)算：
? ? ?Gibbs抽樣與單分量的M-H算法不同之處在于Gibbs抽樣不會(huì)在一些樣本上做停留，即抽樣不會(huì)被拒絕。Gibbs抽樣適用于滿條件分布容易抽樣的情況，而單分量的M-H算法適用于滿條件分布不易抽樣的情況，此時(shí)可選擇容易抽樣的條件分布作為建議分布。
下面以二維正態(tài)分布為例給出多元隨機(jī)變量的Gibbs采樣Python實(shí)現(xiàn)過程。
import math # 導(dǎo)入多元正態(tài)分布函數(shù) from scipy.stats import multivariate_normal # 指定二元正態(tài)分布均值和協(xié)方差矩陣 samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])# 定義給定x的條件下y的條件狀態(tài)轉(zhuǎn)移分布 def p_yx(x, m1, m2, s1, s2):return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2)) # 定義給定y的條件下x的條件狀態(tài)轉(zhuǎn)移分布 def p_xy(y, m1, m2, s1, s2):return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))# 指定相關(guān)參數(shù) N, K= 5000, 20 x_res = [] y_res = [] z_res = [] m1, m2 = 5, -1 s1, s2 = 1, 2 rho, y = 0.5, m2# 遍歷迭代 for i in range(N):for j in range(K):# y給定得到x的采樣x = p_xy(y, m1, m2, s1, s2) # x給定得到y(tǒng)的采樣y = p_yx(x, m1, m2, s1, s2) z = samplesource.pdf([x,y])x_res.append(x)y_res.append(y)z_res.append(z) # 繪圖 num_bins = 50 plt.hist(x_res, num_bins, normed=1, facecolor='green', alpha=0.5,label='x') plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5,label='y') plt.title('Histogram') plt.legend() plt.show();采樣效果如下：
二維效果展示：
MCMC與貝葉斯推斷
? ? ?MCMC對(duì)高效貝葉斯推斷有著重要的作用。假設(shè)有如下貝葉斯后驗(yàn)分布推導(dǎo)：
? ? ?在概率分布多元且形式復(fù)雜的情形下，經(jīng)過貝葉斯先驗(yàn)和似然推導(dǎo)后（即右邊的分式），很難進(jìn)行積分運(yùn)算。具體包括以下三種積分運(yùn)算：規(guī)范化、邊緣化和數(shù)學(xué)期望。
? ? ?首先是上式中的分母，即規(guī)范化計(jì)算：
? ? ?如果是多元隨機(jī)變量或者是包含隱變量，后驗(yàn)分布還需要邊緣化計(jì)算 ：
? ? ?另外如有一個(gè)函數(shù)，可以計(jì)算該函數(shù)關(guān)于后驗(yàn)概率分布的數(shù)學(xué)期望：
? ? ? 當(dāng)觀測(cè)數(shù)據(jù)、先驗(yàn)分布和似然函數(shù)都比較復(fù)雜的時(shí)候，以上三個(gè)積分計(jì)算都會(huì)變得極為困難，這也是早期貝葉斯推斷受到冷落的一個(gè)原因。后來(lái)MCMC方法興起，Bayesian+MCMC的一套做法逐漸流行起來(lái)。
參考資料：
統(tǒng)計(jì)學(xué)習(xí)方法第二版
https://zhuanlan.zhihu.com/p/37121528
往期精彩：
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法20：LDA線性判別分析
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法19：PCA降維
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法18：奇異值分解SVD
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法17：XGBoost
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法16：Adaboost
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法15：GBDT
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法14：Ridge嶺回歸
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法13：Lasso回歸
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法12：貝葉斯網(wǎng)絡(luò)
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法11：樸素貝葉斯
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法10：線性不可分支持向量機(jī)
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法8-9：線性可分支持向量機(jī)和線性支持向量機(jī)
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法7：神經(jīng)網(wǎng)絡(luò)
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法6：感知機(jī)
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法5：決策樹之CART算法
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法4：決策樹之ID3算法
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法3：k近鄰
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法2：邏輯回歸
數(shù)學(xué)推導(dǎo)+純Python實(shí)現(xiàn)機(jī)器學(xué)習(xí)算法1：線性回歸
往期精彩回顧適合初學(xué)者入門人工智能的路線及資料下載機(jī)器學(xué)習(xí)及深度學(xué)習(xí)筆記等資料打印機(jī)器學(xué)習(xí)在線手冊(cè)深度學(xué)習(xí)筆記專輯《統(tǒng)計(jì)學(xué)習(xí)方法》的代碼復(fù)現(xiàn)專輯 AI基礎(chǔ)下載機(jī)器學(xué)習(xí)的數(shù)學(xué)基礎(chǔ)專輯獲取一折本站知識(shí)星球優(yōu)惠券，復(fù)制鏈接直接打開：https://t.zsxq.com/yFQV7am本站qq群1003271085。加入微信群請(qǐng)掃碼進(jìn)群：總結(jié)
以上是生活随笔為你收集整理的【机器学习基础】数学推导+纯Python实现机器学习算法21：马尔可夫链蒙特卡洛...的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： KDD CUP 2020之Debiasi
- 下一篇： 【机器学习基础】数学推导+纯Python