贝叶斯方法学习笔记（二）

原文来自于知乎《你一定从未看过如此通俗易懂的马尔科夫链蒙特卡洛方法解读》。在这里只是将文章主要的核心知识点进行提炼汇总。

链接：你一定从未看过如此通俗易懂的马尔科夫链蒙特卡洛方法解读

MC的基本知识

贝叶斯方法学习笔记（二）
- 一.马尔可夫链及其应用
- 1.马尔可夫链定义
- 2.状态转移矩阵和状态分布
- 3.平稳分布及其存在条件
- 4.马尔可夫链的性质
- 5.可逆马尔可夫链
二.马尔科夫链蒙特卡洛方法(MCMC)
- 1.马尔可夫链蒙特卡洛方法的基本知识
- 2.MCMC的基本步骤
- 3.构造方法
- 4.Metropolis-Hastings算法
- 5.Metropolis-Hastings的代码实现：

一.马尔可夫链及其应用

在蒙特卡洛方法中需要大量采样，构成合适的概率模型，而且需要对采样的随机变量要求其服从某些特定的概率分布，但是有时候满足这种随机分布的随机数是很难产生的，而马尔科夫链蒙特卡洛方法(MCMC)是一个很好的算法。

1.马尔可夫链定义

考虑一个随机变量序列X={X0,X1,X2,...Xt,...}X=\{X_0,X_1,X_2,...X_t,...\}X={X0,X1,X2,...Xt,...},XtX_tXt表示ttt时刻的随机变量。每个XtX_tXt的取值相同，称为状态空间，表示为SSS。既可以连续，又可以离散。当满足以下特性时称为马尔可夫链：

1.当t≥1t \geq1t≥1时，如果XtX_tXt仅仅只依赖X0,X1X_0,X_1X0,X1…,Xt−1X_{t-1}Xt−1，而与{X0,X1，...Xt−2}\{X_0,X_1，...X_{t-2}\}{X0,X1，...Xt−2}无关，即：
P(Xt∣Xt−1,Xt−2,...Xt,...X0)=P(Xt∣Xt−1),t=1,2..P(X_t|X_{t-1},X_{t-2},...X_t,...X_0)=P(X_t|X_{t-1}),t=1,2.. P(Xt∣Xt−1,Xt−2,...Xt,...X0)=P(Xt∣Xt−1),t=1,2..
2.与此同时，在初始时刻的随机变量X0X_0X0遵循P(X0)=π0P(X_0)=\pi_0P(X0)=π0。

称{X0,X1,...Xt,...}\{X_0,X_1,...X_t,...\}{X0,X1,...Xt,...}为马尔可夫链，而P(Xt∣Xt−1)P(X_t|X_{t-1})P(Xt∣Xt−1)称为马尔可夫链的转移概率分布。

当P(Xt∣Xt−1)P(X_t|X_{t-1})P(Xt∣Xt−1)与具体的时间ttt无关时，称为时间齐次的马尔可夫链。

2.状态转移矩阵和状态分布

1.定义Pm×nP_{m\times n}Pm×n为状态转移矩阵，其中pij=P(Xt=i∣Xt−1=j)p_{ij}=P(X_t = i|X_{t-1}=j)pij=P(Xt=i∣Xt−1=j)

2.定义马尔可夫链在ttt时刻的概率分布称为该时刻的状态分布：π(t)=(π0(t),π1(t),...,πn(t))T\pi(t)=(\pi_0(t),\pi_1(t),...,\pi_n(t))^Tπ(t)=(π0(t),π1(t),...,πn(t))T,其中有：
πi(t)=P(Xt=i),i=1,2...\pi_i(t)=P(X_t=i),i=1,2... πi(t)=P(Xt=i),i=1,2...
3.不同时刻的状态分布呈现如下迭代关系：
π(t+1)=Pπ(t)\pi(t+1)=P\pi(t) π(t+1)=Pπ(t)

3.平稳分布及其存在条件

定义当时间ttt足够长时，马尔可夫链的状态分布会收敛到一个常向量上，此时的状态分布称为平稳分布。
给出一个马尔可夫链X={X0,X1,...Xt,...}X=\{X_0,X_1,...X_t,...\}X={X0,X1,...Xt,...}，ttt时刻的状态分布π=(π1,π2,...πn)\pi=(\pi_1,\pi_2,...\pi_n)π=(π1,π2,...πn)是XXX的平稳分布的条件是,π\piπ是方程π=Pπ\pi=P\piπ=Pπ的解，即满足：
xi=∑jpijxj,i=1,2...n,xi,j≥0∑ixi=1x_i=\sum_{j}p_{ij}x_j,i=1,2...n,x_{i,j} \geq 0\\ \sum_ix_i = 1 xi=j∑pijxj,i=1,2...n,xi,j≥0i∑xi=1

例如：当π(0)=(0.5,0.3,0.2)T\pi(0)=(0.5,0.3,0.2)^Tπ(0)=(0.5,0.3,0.2)T，转移矩阵为P=(0.5,0.5,0.250.25,0,0.250.25,0.5,0.5)P=\begin{pmatrix}0.5,0.5,0.25\\0.25,0,0.25\\0.25,0.5,0.5\end{pmatrix}P=⎝⎛0.5,0.5,0.250.25,0,0.250.25,0.5,0.5⎠⎞，求平稳分布：

import numpy as np
import matplotlib.pyplot as plt
start_matrix = np.array([[0.5],[0.3],[0.2]])
transfer_matrix = np.array([[0.5,0.5,0.25],[0.25,0,0.25],[0.25,0.5,0.5]])
init_matrix = start_matrix
next_matrix = []
eps = 1e-5
error = 1
value1_matrix = []
value2_matrix = []
value3_matrix = []
Numbers = 0
while error>eps:next_matrix = np.dot(transfer_matrix,init_matrix)error = np.linalg.norm(next_matrix - init_matrix)init_matrix = next_matrixvalue1_matrix.append(init_matrix[0][0])value2_matrix.append(init_matrix[1][0])value3_matrix.append(init_matrix[2][0])Numbers += 1
print(Numbers)
print(init_matrix)
plt.plot(range(0,Numbers),value1_matrix,'blue',label = 'first')
plt.plot(range(0,Numbers),value2_matrix,'red',label = 'second')
plt.plot(range(0,Numbers),value3_matrix,'yellow',label = 'third')
plt.xlabel("Numbers")
plt.ylabel("Probability")
plt.title("The Markov Chain")
plt.legend()
plt.show()

4.马尔可夫链的性质

1.不可约：从任意状态出发，经过充分长时间后，可以到达任意状态。

给定一个马尔可夫链X={X0,X1,...Xt,...}X=\{X_0,X_1,...X_t,...\}X={X0,X1,...Xt,...},对于∀\forall∀状态i,j∈Si,j\in Si,j∈S,∃t,P(Xt=i∣X0=j)>0\exists t,P(X_t=i|X_0=j)>0∃t,P(Xt=i∣X0=j)>0。

2.非周期：

给定一个马尔可夫链X={X0,X1,...Xt,...}X=\{X_0,X_1,...X_t,...\}X={X0,X1,...Xt,...},对于状态i∈Si\in Si∈S,从t=0t=0t=0,状态iii出发，ttt时刻返回状态的所有时间长{t:P(Xt=i∣X0=i)>0}\{t:P(X_t=i|X_0=i)>0\}{t:P(Xt=i∣X0=i)>0}的最大公约数为1.

3.遍历定理：

不可约且非周期的有限状态马尔可夫链，有唯一平稳分布存在。

5.可逆马尔可夫链

定义：给定马尔可夫链X={X0,X1,...Xt,...}X=\{X_0,X_1,...X_t,...\}X={X0,X1,...Xt,...}，如果有状态分布(π1,π2,..)(\pi_1,\pi_2,..)(π1,π2,..)。∀i,j∈S,∀t\forall i,j\in S,\forall t∀i,j∈S,∀t满足：
P(Xt=i∣Xt−1=j)πj=P(Xt=j∣Xt−1=i)πi,i=1,2,...P(X_t=i|X_{t-1}=j)\pi_j=P(X_t=j|X_{t-1}=i)\pi_i,i=1,2,... P(Xt=i∣Xt−1=j)πj=P(Xt=j∣Xt−1=i)πi,i=1,2,...
称之为XXX的可逆马尔可夫链，上式为细致平衡方程。

定理：满足细致平衡方程的状态分布就是马氏链的平稳分布。

二.马尔科夫链蒙特卡洛方法(MCMC)

1.马尔可夫链蒙特卡洛方法的基本知识

遇到多元变量的随机分布以及复杂的概率密度，需要用马尔可夫链蒙特卡洛方法(MCMC)。

每个时刻在这个马尔可夫链上进行随机游走一次，就可以得到一个样本。根据遍历定理，当时间趋于无穷时，样本的分布趋近平稳分布，样本的函数均值趋近函数的数学期望。

当时间足够长时(t>m)(t>m)(t>m)时,在之后的时间(t<n)(t<n)(t<n)随机游走可以得到样本集合{Xm+1,Xm+2,...Xn}\{X_{m+1},X_{m+2},...X_n\}{Xm+1,Xm+2,...Xn}就是目标概率的抽样结果，而遍历之后的均值就是其数学期望：
Ef=1n−m∑i=m+1nf(Xi)Ef = \frac{1}{n-m}\sum_{i=m+1}^{n}f(X_i) Ef=n−m1i=m+1∑nf(Xi)
关键问题是：

如果给了目标分布p(x)p(x)p(x),相关的马尔可夫链怎么构造？
假设收敛的步数为mmm,迭代mmm步后收敛，总迭代数nnn。m,n=?m,n=?m,n=?

2.MCMC的基本步骤

Step1：构造一个MCMCMC，使得其平稳分布为目标分布p(x)p(x)p(x)。

Step2:从初始分布$ X_0,用构造的,用构造的,用构造的MC$游走，产生样本序列：
{X0,X1,...Xt,...}\{X_0,X_1,...X_t,...\}{X0,X1,...Xt,...}

Step3:由遍历定理，确定m,nm,nm,n，得到平稳样本集合：{Xm+1,Xm+2,...Xt,...Xn}\{X_{m+1},X_{m+2},...X_t,...X_n\}{Xm+1,Xm+2,...Xt,...Xn},求得f(x)f(x)f(x)的均值EfEfEf,由辛钦大数定律：
Ef=1n−m∑i=m+1nf(Xi)Ef = \frac{1}{n-m}\sum_{i=m+1}^{n}f(X_i) Ef=n−m1i=m+1∑nf(Xi)

3.构造方法

构造这样一个MCMCMC的关键点是求出状态转移矩阵PPP

因为满足细致平衡方程的状态分布就是该马尔可夫链的平稳分布Pπ=πP\pi=\piPπ=π，所以找到矩阵PPP即可，即满足：
πiPji=πjPij\pi_iP_{ji}=\pi_jP_{ij} πiPji=πjPij
所以最关键的问题是找到一个匹配上述等式的PPP。

若对于一个任意的马尔可夫状态转移矩阵QijQ_{ij}Qij，构造一个αij\alpha_{ij}αij和αji\alpha_{ji}αji,可以使得：
πiQjiαji=πjQijαij\pi_iQ_{ji}\alpha_{ji}=\pi_jQ_{ij}\alpha_{ij} πiQjiαji=πjQijαij
要使得上式子恒成立，就要取：
αji=πjQijαij=πiQji\alpha_{ji}=\pi_jQ_{ij}\\ \alpha_{ij}=\pi_iQ_{ji} αji=πjQijαij=πiQji
这样可以得到满足细致平衡方程的PPP:
Pji=Qjiαji=πjQijQjiPij=Qijαij=πjQijQjiP_{ji}=Q_{ji}\alpha_{ji}=\pi_jQ_{ij}Q_{ji}\\ P_{ij}=Q_{ij}\alpha_{ij}=\pi_jQ_{ij}Q_{ji} Pji=Qjiαji=πjQijQjiPij=Qijαij=πjQijQji
αij:\alpha_{ij}:αij:接受率，取值在[0,1][0,1][0,1]之间，为一个概率值。

QQQ的平稳分布为建议分布,可以根据上式求出状态转移矩阵PPP

4.Metropolis-Hastings算法

改进的MCMCMCMCMCMC方法：

Step1:Input:Input:Input:状态转移矩阵 QQQ,平稳分布π\piπ,转移次数阈值mmm，需要的样本个数nnn。

Step2:采样一个初始状态X0X_0X0

Step3:forforfor i=0→(n+m)i = 0\rightarrow(n+m)i=0→(n+m)

从MCMCMC (Q)(Q)(Q)中游走一次采样一个x∗x_*x∗

产生随机数u∼U(0,1)u\sim U(0,1)u∼U(0,1)

ififif u<αx∗,xi=min{1,πx∗Qxi,x∗πxiQx∗，xi}u<\alpha_{x_*,x_i}=min\{ 1,\frac{\pi_{x_*}Q_{x_i,x_*}}{\pi_{x_i}Q_{x_*，x_i}}\}u<αx∗,xi=min{1,πxiQx∗，xiπx∗Qxi,x∗}

xi→x∗x_i\rightarrow x_*xi→x∗

xi+1=x∗x_{i+1}=x_*xi+1=x∗

else:else:else:

xi+1=xix_{i+1}=x_ixi+1=xi

output:output:output:样本(xm,xm+1,...xn+m−1)(x_m,x_{m+1},...x_{n+m-1})(xm,xm+1,...xn+m−1)即是平稳分布的样本集

output:output:output:相应的Efp(x)=1n−m∑i=m+1nf(xi)Ef_{p(x)}=\frac{1}{n-m}\sum_{i=m+1}^nf(x_i)Efp(x)=n−m1∑i=m+1nf(xi)

如果选取对称的状态转移矩阵：可以简化αj,i=min{1,πjπi}\alpha_{j,i}=min\{1,\frac{\pi_j}{\pi_i}\}αj,i=min{1,πiπj}

5.Metropolis-Hastings的代码实现：

如果要生成一个a=2.3,b=0.6a=2.3,b=0.6a=2.3,b=0.6的β\betaβ分布。
f(x,a,b)=Γ(a+b)xa−1(1−x)b−1Γ(a)Γ(b)f(x,a,b)=\frac{\Gamma(a+b)x^{a-1}(1-x)^{b-1}}{\Gamma(a)\Gamma(b)} f(x,a,b)=Γ(a)Γ(b)Γ(a+b)xa−1(1−x)b−1
其中Qj,iQ_{j,i}Qj,i是以iii为均值，方差为111的正态分布在jjj处的概率密度的值：

import random
import math
from scipy.stats import beta
from scipy.stats import norm
import matplotlib.pyplot as plta, b = 2, 0.6def norm_dist_prob(theta):y = beta(a, b).pdf(theta)return yT = 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, beta(a, b).pdf(pi),label='Target Distribution', c= 'red')
num_bins = 50
plt.hist(pi, num_bins, density=1, facecolor='green', alpha=0.7,label='Samples Distribution')
plt.legend()
plt.show()

代码来自：你一定从未看过如此通俗易懂的马尔科夫链蒙特卡洛方法解读