Metropolis 采样算法

最新推荐文章于 2025-09-18 12:23:38 发布

原创最新推荐文章于 2025-09-18 12:23:38 发布 · 2.3w 阅读

145 ·

CC 4.0 BY-SA版权

文章标签：

#datamining #python

python 同时被 2 个专栏收录

7 篇文章

订阅专栏

datamining

2 篇文章

订阅专栏

本文介绍了Metropolis 采样算法，这是一种利用马尔科夫链从复杂目标分布中获取样本的方法。首先解释了马尔科夫链的概念、转移概率矩阵和平稳分布，接着详细阐述了Metropolis 采样算法的思路、步骤和实例，展示如何通过该算法对t分布进行采样。文章还探讨了为什么选择对称分布作为基础，并提供了参考文献以深入理解MCMC算法。

部署运行你感兴趣的模型镜像

前言
马尔科夫链
Metropolis 采样算法简介
- Metropolis 采样算法思路
- Metropolis 采样算法步骤
Metropolis 采样算法实例
后记
参考文献

前言

Metropolis 采样算法解决的问题是：从一个复杂的目标分布获取近似的样本。

马尔科夫链

概念描述

$X$ 是一个随机变量，其可能的取值来自于集合 $\Omega=\{s_0,s_1,s_2,\cdots,s_{k-1},s_k\}$ ， $X$ 在离散时刻 t的取值为 $X_t$ 。若 $X$ 随时间变化的转移概率仅仅依赖于其当前时刻的取值 $X_t$ ，即

P (X t + 1 | [X 0, X 1, \dots, X t]) = P (X t + 1 | X t)

$P(X_{t+1}|[X_0,X_1,\cdots,X_t])=P(X_{t+1}|X_t)$
那么随机变量

X X $X$ 随时间变化的过程是一个马尔科夫过程，

X

$X$ 在[0，t]时间内随时间变化生成的序列

(X0,X1,⋯,Xt) ( X 0 , X 1 , ⋯ , X t ) $(X_0,X_1,\cdots,X_t)$ 就是一个马尔科夫链。

转移概率矩阵

设随机变量 $X$ 在任意时刻t+1的取值为 $s_i$ 的概率为 $\pi_i^{t+1}$ ，即 $\pi_i^{t+1}=P(X_{t+1}=s_i)$ ，其中t为一个任意时刻。随机变量从状态 $s_j$ 转移到状态 $s_i$ 的转移概率为 $P_{ji}$ ：

P j i = P (X t + 1 = s i | X t = s j)

$P_{ji}=P(X_{t+1}=s_i|X_t=s_j)$
t 为任意时刻，那么我们可以得到

πt+1i π i t + 1 $\pi_i^{t+1}$ 的计算公式如下：

π t + 1 i = \sum j = 0 k P (X t + 1 = s i | X t = s j) * P (X t = s j)

$\pi_i^{t+1}=\sum_{j=0}^k P(X_{t+1}=s_i|X_t=s_j)*P(X_t=s_j)$
我们注意到在上式中

P(Xt+1=si|Xt=sj)=Pji P ( X t + 1 = s i | X t = s j ) = P j i $P(X_{t+1}=s_i|X_t=s_j)=P_{ji}$ ，

P(Xt=sj)=πtj P ( X t = s j ) = π j t $P(X_t=s_j)=\pi_j^t$ ，因此上式可以写成以下形式：

π t + 1 i = \sum j = 0 k P j i * π t j

$\pi_i^{t+1}=\sum_{j=0}^k P_{ji}*\pi_j^t$
那么我们可以将随机变量

X X $X$ 在t+1时刻取值为

Ω = {s_{0}, s_{1}, s_{2}, \dots, s_{k - 1}, s_{k}}

$\Omega=\{s_0,s_1,s_2,\cdots,s_{k-1},s_k\}$ 的概率写为如下形式：

(π t + 1 0, π t + 1 1, \dots, π t + 1 k) = (π t 0, π t 1, \dots, π t k) * ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ P 00 P 10 ⋮ P k 0 P 01 P 11 ⋮ P k 1 \dots P 0 k \dots P 11 ⋱ ⋮ \dots P k k ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟

$(\pi_0^{t+1},\pi_1^{t+1},\cdots,\pi_k^{t+1})= (\pi_0^t,\pi_1^t,\cdots,\pi_k^t)*\begin{pmatrix}P_{00}&P_{01}&\cdots P_{0k}\\P_{10}&P_{11}&\cdots P_{11}\\\vdots&\vdots&\ddots\vdots\\P_{k0}&P_{k1}&\cdots P_{kk} \end{pmatrix}$
将上式右侧的矩阵记为

Ptran P t r a n $P_{tran}$ ，即

P t r a n = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ P 00 P 10 ⋮ P k 0 P 01 P 11 ⋮ P k 1 \dots P 0 k \dots P 11 ⋱ ⋮ \dots P k k ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟

$P_{tran}=\begin{pmatrix}P_{00}&P_{01}&\cdots P_{0k}\\P_{10}&P_{11}&\cdots P_{11}\\\vdots&\vdots&\ddots\vdots\\P_{k0}&P_{k1}&\cdots P_{kk} \end{pmatrix}$

Ptran P t r a n $P_{tran}$ 就是马尔科夫过程的转移概率矩阵。

平稳分布

如果一个马尔科夫过程的某个状态经过有限次状态转移后又回到自身，那么我们称这个马尔科夫过程具有周期性。如果一个马尔科夫过程存在两个状态，它们之间是互相转移的，那么我们称这个马尔科夫过程可约。若一个马尔科夫过程既没有周期性又不可约，那么这个马尔科夫过程是各态遍历的。各态遍历这个概念我个人的理解是，每个状态都有一定的概率会出现。
各态遍历的马尔科夫过程有一个重要性质：无论初始的状态概率 $(\pi_0^0,\pi_1^0,\cdots,\pi_k^0)$ 取值如何，经过足够多次的状态转移，各态遍历的马尔科夫过程会趋于一个平稳分布 $\pi^s=(\pi_0^s,\pi_1^s,\cdots,\pi_k^s)$ ：

π s = lim t \to \infty π 0 * P t t r a n

$\pi^s=\lim_{t\to\infty}\pi^0*P_{tran}^t$
这时随机变量

X X $X$ 从

Ω = {s_{0}, s_{1}, s_{2}, \dots, s_{k - 1}, s_{k}}

$\Omega=\{s_0,s_1,s_2,\cdots,s_{k-1},s_k\}$ 取值的概率不再变化，而是服从平稳分布

πs π s $\pi^s$ 。
平稳分布满足性质：

πs=πs∗Ptran π s = π s ∗ P t r a n $\pi^s=\pi^s*P_{tran}$ ，即平稳分布经过状态转移后随机变量取值的概率不变，这是当然的。
平稳分布的一个充分条件是：

πi∗Pij=πj∗Pji π i ∗ P i j = π j ∗ P j i $\pi_i*P_{ij}=\pi_j*P_{ji}$ ，对于一个转移概率矩阵

Ptran P t r a n $P_{tran}$ ，若分布

π(X) π ( X ) $\pi(X)$ 满足这个条件，那么它就是

Ptran P t r a n $P_{tran}$ 的平稳分布。这个条件称为 细致平稳条件。

为何要用马尔科夫链

经过上述讨论，有了马尔科夫链的知识基础，我们可以思考：完成上面提出的从一个复杂的目标分布获取近似的样本的任务，马尔科夫链起到了什么作用。我们知道从比较复杂的分布中采样是困难的，如何获得目标分布的样本呢？如果我们先构造出以目标分布为平稳分布的转移概率矩阵 $P_{tran}$ ，然后通过足够多次的状态转移获得了样本序列 $(X_0,X_1,\cdots,X_{n-1},X_n,\cdots)$ ，若在第n次状态转移时马尔科夫过程收敛，这时获得的样本序列 $(X_n,X_{n+1},\cdots)$ 就是从目标分布中采集到的样本，从而间接解决了问题。

Metropolis 采样算法简介

Metropolis 采样算法思路

Metropolis 采样的基本思路是，从一个已知的形式较为简单的分布中采样，并以一定的概率接受这个样本作为目标分布的近似样本。假设我们需要采集一个复杂分布 $\pi(x)$ 的样本，转移概率矩阵 $P$ 以 $\pi(x)$ 为平稳分布。现在我们有一个形式较简单的对称分布 $\theta(x)$ （下文中会说明为什么一定要选择对称分布），转移概率矩阵 $Q$ 以其为平稳分布，从 $\theta(x|x_{n-1})$ 中采样得到一个候选样本 $\hat x$ ，以一定的概率选择接受或拒绝这个候选样本。一般这个接受的概率定为 $\alpha=min(1,\frac{\pi(\hat x)}{\pi(x_{n-1})})$ ，当随机生成的一个概率值小于这个指定的概率时，接受这个候选样本，即 $x_n=\hat x$ ，否则拒绝这个候选样本并令 $x_n=x_{n-1}$ 。
从分布 $\theta(x|x_{n-1})$ 中采出的样本为什么可以作为 $\pi(x)$ 的样本的近似呢？下面给出一点理论上的支撑。 $\theta(x)$ 的转移概率矩阵为 $Q$ ，随机变量从状态 $s_i$ 转移到 $s_j$ 的转移概率为 $Q_{ij}$ ，随机变量从状态 $s_i$ 转移到 $s_j$ 的接受概率为 $\alpha_{ij}$ ，我们借助于 $Q_{ij}和\alpha_{ij}$ 得到的转移概率 $P_{ij}和Q_{ij}、\alpha_{ij}$ 有如下关系：

P i j = α i j * Q i j

$P_{ij}=\alpha_{ij}*Q_{ij}$
我们只需要保证转移概率矩阵

P P $P$ 满足细致平稳条件

π (i) P_{i j} = π (j) P_{j i}

$\pi(i)P_{ij}=\pi(j)P_{ji}$ ，那么

P P $P$ 的平稳分布就是目标分布

π (x)

$\pi(x)$ ，马尔科夫链收敛后的样本也就是目标分布

π(x) π ( x ) $\pi(x)$ 的样本。下面证明转移概率矩阵

P P $P$ 满足细致平稳条件

π (i) P_{i j} = π (j) P_{j i}

$\pi(i)P_{ij}=\pi(j)P_{ji}$ ：

π (i) P i j = π (i) * α i j * Q i j = π (i) * m i n (1, π ( j ) π ( i )) * Q i j

$\pi(i)P_{ij}=\pi(i)*\alpha_{ij}*Q_{ij}=\pi(i)*min(1,\frac{\pi(j)}{\pi(i)})*Q_{ij}$

π (i) P i j = m i n (π (i) * Q i j, π (j) * Q i j)

$\pi(i)P_{ij}=min(\pi(i)*Q_{ij},\pi(j)*Q_{ij})$
在上文中我们提到

θ(x) θ ( x ) $\theta(x)$ 为对称分布，因此

Qij=Qji Q i j = Q j i $Q_{ij}=Q_{ji}$ 。上面的式子就变成：

π (i) P i j = m i n (π (i) * Q j i, π (j) * Q j i) = π (j) * m i n (1, π ( i ) π ( j )) * Q j i

$\pi(i)P_{ij}=min(\pi(i)*Q_{ji},\pi(j)*Q_{ji})=\pi(j)*min(1,\frac{\pi(i)}{\pi(j)})*Q_{ji}$

π (i) P i j = π (j) * α j i * Q j i = π (j) * P j i

$\pi(i)P_{ij}=\pi(j)*\alpha_{ji}*Q_{ji}=\pi(j)*P_{ji}$
因此转移概率矩阵

P P $P$ 的平稳分布为

π (x)

$\pi(x)$ ，以上述方式获得的马尔科夫链收敛后的样本序列为目标分布

π(x) π ( x ) $\pi(x)$ 的近似样本。

Metropolis 采样算法步骤

Metropolis 采样算法
step1. 初始化：t=0，随机生成一个 $x_0$ 赋值给当前的 $x_t$ ，迭代终止条件为 t=T
step2. 令t=t+1，从条件概率分布 $\theta(x\|x_{t-1})$ 中生成候选样本 $\hat x$
step3. 计算接受概率 $\alpha$ ， $\alpha=min（1，\frac{\pi(\hat x)}{\pi(x_{t-1})}）$
step4. 从均匀分布中生成一个随机数 $\alpha_t$
step5. 若 $\alpha_t\le\alpha$ ，则接受候选样本， $x_t=\hat x$ ，否则，拒绝候选样本并令 $x_t=x_{t-1}$
step6. 若t=T，停止迭代，否则回到第2步继续迭代

算法停止迭代后，生成了样本序列 $(x_0,x_1,\cdots,x_{T-1},x_T)$ ，根据需要截取尾部的n个样本 $(\underbrace{x_{T-n+1},x_{T-n+2},\cdots,x_{T-1},x_T}_{n})$ ，这n个样本可以近似地认为是从目标分布 $\pi(x)$ 中采样得到的。

Metropolis 采样算法实例

下面给出利用Metropolis 采样算法对t分布进行采样的一个例子。t分布的概率密度函数为：

f (x, n) = Γ ( n + 1 2 ) n π - - - \sqrt Γ ( n 2 ) (1 + x 2 n) - n + 1 2

$f(x,n)=\frac{\Gamma(\frac{n+1}{2})}{\sqrt{n\pi}\Gamma(\frac{n}{2})}(1+\frac{x^2}{n})^{-\frac{n+1}{2}}$
其中n是t分布的自由度。下图是自由度为3的t分布的概率密度函数分布图。

下面的例子中自由度取3，样本序列的初始值从均匀分布中随机选取，从正态分布中采集候选样本，迭代次数为10000次，取样本序列的第9001到第10000个样本作为从t分布中采集的样本的近似，并画出选取的1000个样本的分布直方图，与标准分布进行比较。

import numpy as np
from scipy.stats import uniform
from scipy.stats import norm
from scipy.stats import t
import matplotlib.pyplot as plt

# 自由度为3的t分布
def t_distribution(x):    # n=3
    p = 2/(np.sqrt(3)*np.pi*np.square(1+np.square(x)/3))
    return p

T = 10000   # 迭代次数
sigma = 1.  # 正态分布标准差
sample_x = np.zeros(T+1)
sample_x[0] = uniform.rvs(size=1)[0]   # 初始化马尔科夫链初值
for i in range (1, T+1):
    hat_x = norm.rvs(loc = sample_x[i-1], scale=sigma, size=1, random_state=None)   # 从正态分布中生成候选值
    alpha = min(1, t_distribution(hat_x[0])/t_distribution(sample_x[i-1]))  # 计算接受概率
    alpha_t = uniform.rvs(size=1)[0]  # 生成接受概率判定值
    if alpha_t <= alpha :      # 若判定值小于接受概率则接受候选值，否则拒绝候选值
        sample_x[i] = hat_x[0]
    else:
        sample_x[i] = sample_x[i-1]

fig, ax = plt.subplots(1, 1)
df = 3   # t分布的自由度为3
mean, var, skew, kurt = t.stats(df, moments='mvsk')
x = np.linspace(t.ppf(0.01, df), t.ppf(0.99, df), 100)
p1 = ax.plot(x, t.pdf(x, df), 'k-', lw=5, alpha=0.6, label='t pdf')     # pdf: Probability density function；画自由度为3的标准t分布曲线
p2 = plt.hist(sample_x[9001:], 100, normed=True, alpha=0.5, histtype='stepfilled', facecolor='red', label='sample_t')   # 画生成的马尔科夫链的标准化柱状图
plt.legend()
plt.show()