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

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

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

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

一.马尔可夫链及其应用

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

1.马尔可夫链定义

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

​ 1.当$t\ge 1$时，如果$X_t$仅仅只依赖$X_0,X_1$…,$X_{t-1}$，而与 { X 0 , X 1 ， . . . X t − 2 } \{X_0,X_1，...X_{t-2}\} {X0​,X1​，...Xt−2​}无关，即：
$$P(X_t|X_{t-1},X_{t-2},...X_t,...X_0)=P(X_t|X_{t-1}),t=1,2..$$
​ 2.与此同时，在初始时刻的随机变量$X_0$遵循$P(X_0)={\pi }_0$。

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

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

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

​ 1.定义$P_{m\times n}$为状态转移矩阵，其中$p_{ij}=P(X_t=i|X_{t-1}=j)$

​ 2.定义马尔可夫链在$t$时刻的概率分布称为该时刻的状态分布：$\pi (t)=({\pi }_0(t),{\pi }_1(t),...,{\pi }_n(t){)}^T$,其中有：
$${\pi }_i(t)=P(X_t=i),i=1,2...$$
​ 3.不同时刻的状态分布呈现如下迭代关系：
$$\pi (t+1)=P\pi (t)$$

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

1. 定义当时间$t$足够长时，马尔可夫链的状态分布会收敛到一个常向量上，此时的状态分布称为平稳分布。

2. 给出一个马尔可夫链 X = { X 0 , X 1 , . . . X t , . . . } X=\{X_0,X_1,...X_t,...\} X={X0​,X1​,...Xt​,...}，$t$时刻的状态分布$\pi =({\pi }_1,{\pi }_2,...{\pi }_n)$是$X$的平稳分布的条件是,$\pi$是方程$\pi =P\pi$的解，即满足：
   $$\begin{gathered} x_i=\sum _j{p}_{ij}{x}_j,i=1,2...n,x_{i,j}\ge 0 \\ \sum _i{x}_i=1 \end{gathered}$$

3. 例如：当$\pi (0)=(0.5,0.3,0.2{)}^T$，转移矩阵为$P=\left(\begin{array}{c}0.5,0.5,0.25\\ 0.25,0,0.25\\ 0.25,0.5,0.5\end{array}\right)$，求平稳分布：

   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_matrix
       value1_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 = { X 0 , X 1 , . . . X t , . . . } X=\{X_0,X_1,...X_t,...\} X={X0​,X1​,...Xt​,...},对于$\forall$状态$i,j\in S$,$\exists t,P(X_t=i|X_0=j)>0$。

​ 2.非周期：

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

​ 3.遍历定理：

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

5.可逆马尔可夫链

定义：给定马尔可夫链 X = { X 0 , X 1 , . . . X t , . . . } X=\{X_0,X_1,...X_t,...\} X={X0​,X1​,...Xt​,...}，如果有状态分布$({\pi }_1,{\pi }_2,..)$。$\forall i,j\in S,\forall t$满足：
$$P(X_t=i|X_{t-1}=j){\pi }_j=P(X_t=j|X_{t-1}=i){\pi }_i,i=1,2,...$$
​ 称之为$X$的可逆马尔可夫链，上式为细致平衡方程。

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

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

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

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

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

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

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

2.MCMC的基本步骤

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

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

​ Step3:由遍历定理，确定$m,n$，得到平稳样本集合： { X m + 1 , X m + 2 , . . . X t , . . . X n } \{X_{m+1},X_{m+2},...X_t,...X_n\} {Xm+1​,Xm+2​,...Xt​,...Xn​},求得$f(x)$的均值$Ef$,由辛钦大数定律：
$$Ef=\frac{1}{n-m}\sum _{i=m+1}^{n}f(X_i)$$

3.构造方法

​ 构造这样一个$MC$的关键点是求出状态转移矩阵$P$

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

​ 若对于一个任意的马尔可夫状态转移矩阵$Q_{ij}$，构造一个${\alpha }_{ij}$和${\alpha }_{ji}$,可以使得：
$${\pi }_i{Q}_{ji}{\alpha }_{ji}={\pi }_j{Q}_{ij}{\alpha }_{ij}$$
​ 要使得上式子恒成立，就要取：
$$\begin{gathered} {\alpha }_{ji}={\pi }_j{Q}_{ij} \\ {\alpha }_{ij}={\pi }_i{Q}_{ji} \end{gathered}$$
​ 这样可以得到满足细致平衡方程的$P$:
$$\begin{gathered} P_{ji}=Q_{ji}{\alpha }_{ji}={\pi }_j{Q}_{ij}{Q}_{ji} \\ {P}_{ij}=Q_{ij}{\alpha }_{ij}={\pi }_j{Q}_{ij}{Q}_{ji} \end{gathered}$$
​ ${\alpha }_{ij}:$接受率，取值在$[0,1]$之间，为一个概率值。

​ $Q$的平稳分布为建议分布,可以根据上式求出状态转移矩阵$P$

4.Metropolis-Hastings算法

改进的$MCMC$方法：

Step1:$Input:$状态转移矩阵 $Q$,平稳分布$\pi$,转移次数阈值$m$，需要的样本个数$n$。

Step2:采样一个初始状态$X_0$

Step3:$for$ $i=0\to (n+m)$

​ 从$MC$ $(Q)$中游走一次采样一个$x_{\ast }$

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

​ $if$ u < α x ∗ , x i = m i n { 1 , π x ∗ Q x i , x ∗ π x i Q x ∗ ， x i } 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,πxi​​Qx∗​，xi​​πx∗​​Qxi​,x∗​​​}

​ $x_i\to x_{\ast }$

​ $x_{i+1}=x_{\ast }$

​ $else:$

​ $x_{i+1}=x_i$

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

$output:$相应的$E{f}_{p(x)}=\frac{1}{n-m}{\sum }_{i=m+1}^{n}f(x_i)$

如果选取对称的状态转移矩阵：可以简化 α j , i = m i n { 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.6$的$\beta$分布。
$$f(x,a,b)=\frac{\Gamma (a+b)x^{a-1}(1-x{)}^{b-1}}{\Gamma (a)\Gamma (b)}$$
其中$Q_{j,i}$是以$i$为均值，方差为$1$的正态分布在$j$处的概率密度的值：

import random
import math
from scipy.stats import beta
from scipy.stats import norm
import matplotlib.pyplot as plt

a, b = 2, 0.6

def norm_dist_prob(theta):
    y = beta(a, b).pdf(theta)
    return y

T = 5000
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, (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()

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