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

本文是李航老師的《統(tǒng)計學(xué)習(xí)方法》一書的代碼復(fù)現(xiàn)。作者：黃海廣

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

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

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

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

計算樣本均值，作為數(shù)學(xué)期望的估計值。

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

通常考慮時間齊次馬爾可夫鏈。有離散狀態(tài)馬爾可夫鏈和連續(xù)狀態(tài)馬爾可夫鏈，分別由概率轉(zhuǎn)移矩陣和概率轉(zhuǎn)移核定義。

滿足或的狀態(tài)分布稱為馬爾可夫鏈的平穩(wěn)分布。

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

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

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

（1）在隨機變量的狀態(tài)空間上構(gòu)造一個滿足遍歷定理條件的馬爾可夫鏈，其平穩(wěn)分布為目標分布；

（2）由狀態(tài)空間的某一點出發(fā)，用所構(gòu)造的馬爾可夫鏈進行隨機游走，產(chǎn)生樣本序列；

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

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

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

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

---

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

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

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

之后使用該平穩(wěn)分布的樣本進行近似數(shù)值計算。

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

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

馬爾可夫鏈蒙特卡羅法被應(yīng)用于概率分布的估計、定積分的近似計算、最優(yōu)化問題的近似求解等問題，特別是被應(yīng)用于統(tǒng)計學(xué)習(xí)中概率模型的學(xué)習(xí)

與推理，是重要的統(tǒng)計學(xué)習(xí)計算方法。

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

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

19.1.2 數(shù)學(xué)期望估計

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

假設(shè)有隨機變量， 取值 , 其概率密度函數(shù)為 , 為定義在 上的函數(shù)， 目標是求函數(shù) 關(guān)于密度函數(shù) 的數(shù)學(xué)期望 。

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

數(shù)的樣本均值

作為數(shù)學(xué)期望近似值。

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

這樣就得到了數(shù)學(xué)期望的近似計算方法：

馬爾可夫鏈

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

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

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

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

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

，， 這一性質(zhì)稱為馬爾可夫性，即

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

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

平穩(wěn)分布

設(shè)有馬爾可夫鏈，其狀態(tài)空間為 ,轉(zhuǎn)移概率矩陣為 ， 如果存在狀態(tài)空間 上的一個分布

使得

則稱丌為馬爾可夫褳的平穩(wěn)分布。

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

引理19．1

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

吉布斯采樣

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

輸出： 的隨機樣本 ，函數(shù)樣本均值 ;

參數(shù)： 收斂步數(shù), 迭代步數(shù) .

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

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

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

• ...

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

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

得到第 次迭代值 .

得到樣本集合

{}

計算

---

網(wǎng)絡(luò)資源：

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()

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

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

M-H采樣python實現(xiàn)

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

假設(shè)目標平穩(wěn)分布是一個均值3，標準差2的正態(tài)分布，而選擇的馬爾可夫鏈狀態(tài)轉(zhuǎn)移矩陣 的條件轉(zhuǎn)移概率是以 為均值,方差1的正態(tài)分布在位置 的值。

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) #狀態(tài)轉(zhuǎn)移進行隨機抽樣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實現(xiàn)

假設(shè)我們要采樣的是一個二維正態(tài)分布 ，其中： ,

;

而采樣過程中的需要的狀態(tài)轉(zhuǎn)移條件分布為：

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(tǒng)的采樣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] 《統(tǒng)計學(xué)習(xí)方法》: https://baike.baidu.com/item/統(tǒng)計學(xué)習(xí)方法/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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。