第19章 馬爾可夫鏈蒙特卡羅法

本文是李航老師的《統計學習方法》一書的代碼復現。作者：黃海廣

備注：代碼都可以在github中下載。我將陸續將代碼發布在公眾號“機器學習初學者”，可以在這個專輯在線閱讀。

蒙特卡羅法是通過基于概率模型的抽樣進行數值近似計算的方法，蒙特卡羅法可以用于概率分布的抽樣、概率分布數學期望的估計、定積分的近似計算。

隨機抽樣是蒙特卡羅法的一種應用，有直接抽樣法、接受拒絕抽樣法等。接受拒絕法的基本想法是，找一個容易抽樣的建議分布，其密度函數的數倍大于等于想要抽樣的概率分布的密度函數。按照建議分布隨機抽樣得到樣本，再按要抽樣的概率分布與建議分布的倍數的比例隨機決定接受或拒絕該樣本，循環執行以上過程。

馬爾可夫鏈蒙特卡羅法數學期望估計是蒙特卡羅法的另一種應用，按照概率分布抽取隨機變量的個獨立樣本，根據大數定律可知，當樣本容量增大時，函數的樣本均值以概率1收斂于函數的數學期望

計算樣本均值，作為數學期望的估計值。

馬爾可夫鏈是具有馬爾可夫性的隨機過程

通常考慮時間齊次馬爾可夫鏈。有離散狀態馬爾可夫鏈和連續狀態馬爾可夫鏈，分別由概率轉移矩陣和概率轉移核定義。

滿足或的狀態分布稱為馬爾可夫鏈的平穩分布。

馬爾可夫鏈有不可約性、非周期性、正常返等性質。一個馬爾可夫鏈若是不可約、非周期、正常返的，則該馬爾可夫鏈滿足遍歷定理。當時間趨于無窮時，馬爾可夫鏈的狀態分布趨近于平穩分布，函數的樣本平均依概率收斂于該函數的數學期望。

可逆馬爾可夫鏈是滿足遍歷定理的充分條件。

馬爾可夫鏈蒙特卡羅法是以馬爾可夫鏈為概率模型的蒙特卡羅積分方法，其基本想法如下：

（1）在隨機變量的狀態空間上構造一個滿足遍歷定理條件的馬爾可夫鏈，其平穩分布為目標分布；

（2）由狀態空間的某一點出發，用所構造的馬爾可夫鏈進行隨機游走，產生樣本序列；

（3）應用馬爾可夫鏈遍歷定理，確定正整數和$n(m<n)$，得到樣本集合$\{ 1="" 2="" x="" _="" {="" m="" +="" }="" ,="" \cdots="" n="" \}$，進行函數$f(x)$的均值（遍歷均值）估計：<="" p="">

Metropolis-Hastings算法是最基本的馬爾可夫鏈蒙特卡羅法。假設目標是對概率分布進行抽樣，構造建議分布，定義接受分布進行隨機游走，假設當前處于狀態，按照建議分布機抽樣，按照概率接受抽樣，轉移到狀態 ，按照概率拒絕抽樣，停留在狀態，持續以上操作，得到一系列樣本。這樣的隨機游走是根據轉移核為的可逆馬爾可夫鏈（滿足遍歷定理條件）進行的，其平穩分布就是要抽樣的目標分布。

吉布斯抽樣（Gibbs sampling）用于多元聯合分布的抽樣和估計吉布斯抽樣是單分量 Metropolis-Hastings算法的特殊情況。這時建議分布為滿條件概率分布

吉布斯抽樣的基本做法是，從聯合分布定義滿條件概率分布，依次從滿條件概率分布進行抽樣，得到聯合分布的隨機樣本。假設多元聯合概率分布為，吉布斯抽樣從一個初始樣本出發，不斷進行迭代，每一次迭代得到聯合分布的一個樣本，在第次迭代中，依次對第個變量按照滿條件概率分布隨機抽樣，得到最終得到樣本序列。

---

蒙特卡洛法（Monte Carlo method) ， 也稱為統計模擬方法 (statistical simulation method) ， 是通過從概率模型的隨機抽樣進行近似數值計

算的方法。 馬爾可夫鏈陟特卡羅法 (Markov Chain Monte Carlo, MCMC)， 則是以馬爾可夫鏈 (Markov chain）為概率模型的蒙特卡洛法。

馬爾可夫鏈蒙特卡羅法構建一個馬爾可夫鏈，使其平穩分布就是要進行抽樣的分布， 首先基于該馬爾可夫鏈進行隨機游走， 產生樣本的序列，

之后使用該平穩分布的樣本進行近似數值計算。

Metropolis-Hastings算法是最基本的馬爾可夫鏈蒙特卡羅法，Metropolis等人在 1953年提出原始的算法，Hastings在1970年對之加以推廣，

形成了現在的形式。吉布斯抽樣(Gibbs sampling)是更簡單、使用更廣泛的馬爾可夫鏈蒙特卡羅法，1984 年由S. Geman和D. Geman提出。

馬爾可夫鏈蒙特卡羅法被應用于概率分布的估計、定積分的近似計算、最優化問題的近似求解等問題，特別是被應用于統計學習中概率模型的學習

與推理，是重要的統計學習計算方法。

一般的蒙特卡羅法有直接抽樣法、接受-拒絕抽樣法、 重要性抽樣法等。

接受-拒絕抽樣法、重要性抽樣法適合于概率密度函數復雜 （如密度函數含有多個變量，各變量相互不獨立，密度函數形式復雜），不能直接抽樣的情況。

19.1.2 數學期望估計

一艤的蒙特卡羅法， 如直接抽樣法、接受·拒絕抽樣法、重要性抽樣法， 也可以用于數學期望估計 （estimation Of mathematical expectation)。

假設有隨機變量， 取值 , 其概率密度函數為 , 為定義在 上的函數， 目標是求函數 關于密度函數 的數學期望 。

針對這個問題，蒙特卡羅法按照概率分布 獨立地抽取 個樣本，比如用以上的抽樣方法，之后計算函

數的樣本均值

作為數學期望近似值。

根據大數定律可知， 當樣本容量增大時， 樣本均值以概率1收斂于數學期望：

這樣就得到了數學期望的近似計算方法：

馬爾可夫鏈

考慮一個隨機變量的序列 這里 ，表示時刻 的隨機變量， .

每個隨機變量 的取值集合相同， 稱為狀態空間， 表示為. ?隨機變量可以是離散的， 也可以是連續的。

以上隨機變量的序列構成隨機過程（stochastic process)。

假設在時刻 的隨機變量 遵循概率分布 ，稱為初始狀態分布。在某個時刻 的隨機變量 與前

一個時刻的隨機變量 之間有條件分布 如果 只依賴于 , 而不依賴于過去的隨機變量

，， 這一性質稱為馬爾可夫性，即

具有馬爾可夫性的隨機序列稱為馬爾可夫鏈， 或馬爾可夫過程（Markov process)。 條件概率分布

稱為馬爾可夫鏈的轉移概率分布。 轉移概率分布決定了馬爾可夫褳的特性。

平穩分布

設有馬爾可夫鏈，其狀態空間為 ,轉移概率矩陣為 ， 如果存在狀態空間 上的一個分布

使得

則稱丌為馬爾可夫褳的平穩分布。

直觀上，如果馬爾可夫鏈的平穩分布存在，那么以該平穩分布作為初始分布，面向未來進行隨機狀態轉移，之后任何一個時刻的狀態分布都是該平穩分布。

引理19．1

給定一個馬爾可夫鏈, 狀態空間為, 移概率矩陣為， 則分布 為 的平穩分布的充要條件是是下列方程組的解：

吉布斯采樣

輸入： 目標概率分布的密度函數, 函數;

輸出： 的隨機樣本 ，函數樣本均值 ;

參數： 收斂步數, 迭代步數 .

初始化。給出初始樣本 ($x^{0}{1}, x^{0}{2},..., x^{0}_{k}^{T}$.

對循環執行
設第次迭代結束前的樣本為($x^{i-1}{1}, x^{i-1}{2},..., x^{i-1}_{k}^{T}，則第i$次迭代進行如下幾步操作：

• (1)由滿條件分布 抽取

• ...

• (j)由滿條件分布 抽取

• (k)由滿條件分布 抽取

得到第 次迭代值 .

得到樣本集合

{}

計算

---

網絡資源：

LDA-math-MCMC 和 Gibbs Sampling： https://cosx.org/2013/01/lda-math-mcmc-and-gibbs-sampling

MCMC蒙特卡羅方法： https://www.cnblogs.com/pinard/p/6625739.html

import random import math import matplotlib.pyplot as plt import seaborn as sns import numpy as nptransfer_matrix = np.array([[0.6, 0.2, 0.2], [0.3, 0.4, 0.3], [0, 0.3, 0.7]],dtype='float32') start_matrix = np.array([[0.5, 0.3, 0.2]], dtype='float32')value1 = [] value2 = [] value3 = [] for i in range(30):start_matrix = np.dot(start_matrix, transfer_matrix)value1.append(start_matrix[0][0])value2.append(start_matrix[0][1])value3.append(start_matrix[0][2]) print(start_matrix) [[0.23076935 0.30769244 0.46153864]] #進行可視化 x = np.arange(30) plt.plot(x,value1,label='cheerful') plt.plot(x,value2,label='so-so') plt.plot(x,value3,label='sad') plt.legend() plt.show()

可以發現，從10輪左右開始，我們的狀態概率分布就不變了，一直保持在
[0.23076934,0.30769244,0.4615386]

參考：https://zhuanlan.zhihu.com/p/37121528

M-H采樣python實現

https://zhuanlan.zhihu.com/p/37121528

假設目標平穩分布是一個均值3，標準差2的正態分布，而選擇的馬爾可夫鏈狀態轉移矩陣 的條件轉移概率是以 為均值,方差1的正態分布在位置 的值。

from scipy.stats import norm def norm_dist_prob(theta):y = norm.pdf(theta, loc=3, scale=2)return y T = 5000 pi = [0 for i in range(T)] sigma = 1 t = 0 while t < T - 1:t = t + 1pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1,random_state=None) #狀態轉移進行隨機抽樣alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1]))) #alpha值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 = 50 plt.hist(pi,num_bins,density=1,facecolor='red',alpha=0.7,label='Samples Distribution') plt.legend() plt.show()

二維Gibbs采樣實例python實現

假設我們要采樣的是一個二維正態分布 ，其中： ,

;

而采樣過程中的需要的狀態轉移條件分布為：

from mpl_toolkits.mplot3d import Axes3D from scipy.stats import multivariate_normalsamplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])def p_ygivenx(x, m1, m2, s1, s2):return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))def p_xgiveny(y, m1, m2, s1, s2):return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))N = 5000 K = 20 x_res = [] y_res = [] z_res = [] m1 = 5 m2 = -1 s1 = 1 s2 = 2rho = 0.5 y = m2for i in range(N):for j in range(K):x = p_xgiveny(y, m1, m2, s1, s2) #y給定得到x的采樣y = p_ygivenx(x, m1, m2, s1, s2) #x給定得到y的采樣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,density=1, facecolor='green', alpha=0.5,label='x') plt.hist(y_res, num_bins, density=1, facecolor='red', alpha=0.5,label='y') plt.title('Histogram') plt.legend() plt.show()

---

本章代碼來源：https://github.com/hktxt/Learn-Statistical-Learning-Method

下載地址

https://github.com/fengdu78/lihang-code

參考資料：

[1] 《統計學習方法》: https://baike.baidu.com/item/統計學習方法/10430179

[2] 黃海廣: https://github.com/fengdu78

[3] ?github: https://github.com/fengdu78/lihang-code

[4] ?wzyonggege: https://github.com/wzyonggege/statistical-learning-method

[5] ?WenDesi: https://github.com/WenDesi/lihang_book_algorithm

[6] ?火燙火燙的: https://blog.csdn.net/tudaodiaozhale

[7] ?hktxt: https://github.com/hktxt/Learn-Statistical-Learning-Method

總結

以上是生活随笔為你收集整理的复现经典：《统计学习方法》第19章 马尔可夫链蒙特卡罗法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。