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