【機器學習基礎】數學推導+純Python實現機器學習算法21:馬爾可夫鏈蒙特卡洛

Python機器學習算法實現

Author：louwill

Machine Learning Lab

     

     蒙特卡洛（Monte Carlo，MC）方法作為一種統計模擬和近似計算方法，是一種通過對概率模型隨機抽樣進行近似數值計算的方法。馬爾可夫鏈（Markov Chain，MC）則是一種具備馬爾可夫性的隨機序列。將二者結合起來便有了馬爾可夫鏈蒙特卡洛方法（Markov Chain Monte Carlo，MCMC），即是以構造馬爾可夫鏈為概率模型的蒙特卡洛方法。     限於篇幅，本文不對MCMC的前置知識進行詳細介紹，關於蒙特卡洛方法和馬爾可夫鏈的大量內容，請各位讀者自行查閱相關材料進行學習。本文聚焦於MCMC方法本身原理和常用實現方法。MCMC簡介     一般來說，對目標概率模型進行隨機抽樣能夠幫助我們得到該分布的近似數值解。但如果隨機變量的多元的，或者所要抽樣的概率密度函數形式是複雜的非標準式時，直接應用蒙特卡洛方法就會很困難。     MCMC方法的基本思路是：在隨機變量

from scipy.stats import normimport randomimport matplotlib.pyplot as plt
def smooth_dist(theta):    y = norm.pdf(theta, loc=3, scale=2)    return y
T = 10000pi = [0 for i in range(T)]sigma = 1t = 0while 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])))       u = 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 = 100plt.hist(pi,          num_bins,          normed=1,          facecolor='red',          alpha=0.6,          label='Samples Distribution')plt.legend()plt.show();
     相較於M-H抽樣，Gibbs抽樣是目前更常用的MCMC抽樣算法，Gibbs抽樣可以視為特殊的M-H抽樣方法。Gibbs抽樣適用於多元隨機變量聯合分布的抽樣和估計，其基本思路是從聯合概率分布種定義滿條件概率分布，依次對滿條件概率分布進行抽樣得到目標樣本序列。     這裡簡單提一下滿條件分布。假設MCMC的目標分布為多元聯合概率分布
import mathfrom scipy.stats import multivariate_normalsamplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])
def p_yx(x, m1, m2, s1, s2):    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))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, 20x_res = []y_res = []z_res = []m1, m2 = 5, -1s1, s2 = 1, 2rho, y = 0.5, m2
for i in range(N):    for j in range(K):                x = p_xy(y, m1, m2, s1, s2)                  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 = 50plt.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();
二維效果展示：
from mpl_toolkits.mplot3d import Axes3Dfig = plt.figure()ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20)ax.scatter(x_res, y_res, z_res, marker='o')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：線性回歸
獲取一折本站知識星球優惠券，複製連結直接打開：
https://t.zsxq.com/yFQV7am
本站qq群1003271085。
加入微信群請掃碼進群：