29、马尔可夫链蒙特卡罗法

healed萌

已于 2023-11-27 16:47:29 修改

阅读量1.1k

点赞数 23

分类专栏：机器学习文章标签：机器学习人工智能

于 2023-11-27 12:30:10 首次发布

本文链接：https://blog.csdn.net/m0_56642803/article/details/134586696

版权

机器学习专栏收录该内容

30 篇文章 0 订阅

订阅专栏

马尔可夫链蒙特卡罗法

蒙特卡罗法（Monte Carlo method），也称为统计模拟方法（statistical simulation method），是通过从概率模型的随机抽样进行 近似数值计算 的方法
马尔可夫链蒙特卡罗法（Markov Chain Monte Carlo，MCMC），则是以马尔可夫链（Markov chain）为概率模型的蒙特卡罗法
马尔可夫链蒙特卡罗法构建一个马尔可夫链，使其平稳分布就是要进行抽样的分布，首先基于该马尔可夫链进行随机游走，产生样本的序列，之后使用该平稳分布的样本进行近似数值计算
马尔可夫链蒙特卡罗法被应用于概率分布的估计、定积分的近似计算、最优化问题的近似求解等问题，特别是被应用于统计学习中概率模型的学习与推理，是重要的统计学习计算方法

1 蒙特卡罗法

1.1 随机抽样

统计学和机器学习的目的是基于数据对概率分布的特征进行推断。蒙特卡罗法要解决的问题是，假设概率分布的定义已知，通过抽样获得概率分布的随机样本，并通过得到的随机样本对概率分布的特征进行分析。

核心是随机抽样。一般的蒙特卡罗法有直接抽样法、接受-拒绝抽样法、重要性抽样法。接受-拒绝抽样法、重要性抽样法适合于概率密度函数复杂（如密度函数含有多个变量，各变量相互不独立，密度函数形式复杂），不能直接抽样的情况。

接受-拒绝抽样法：假设有随机变量 $x$ ，取值 $x \in X$ ，其概率密度函数为 $p (x)$ ；目标是得到该概率分布的随机样本，以对这个概率分布进行分析。基本想法是假设 $p (x)$ 不可以直接抽样，找一个可以直接抽样的分布，称为建议分布；假设 $q (x)$ 是建议分布的概率密度函数，并且有 $q (x)$ 的 $c$ 倍一定大于等于 $p (x)$ ，其中 $c ＞ 0$ ；按照 $q (x)$ 进行抽样，假设得到结果是 $x^*$ ，再按照 $p(x^*)/(cq(x^*))$ 的比例随机决定是否接受 $x^*$ ；直观上，落到 $p(x^*)$ 范围内的就接受，落到 $p(x^*)$ 范围外就拒绝。接受-拒绝法实际上是按照 $p (x)$ 的涵盖面积（或涵盖体积）占 $c q (x)$ 的涵盖面积（或涵盖体积）的比例进行抽样。优点是容易实现，缺点是效率不高。

在这里插入图片描述

接受-拒绝法具体算法如下：

输入：抽样的目标概率分布的概率密度函数 $p (x)$ ;
输出：概率分布的随机样本 $x_1, x_2,… ,x_n$ 。
参数：样本数 $n$
(1) 选择概率密度函数为 $q (x)$ 的概率分布，作为建议分布，使其对任一 $x$ 满足 $\geq p(x)$ ，其中 $c > 0$ 。
(2) 按照建议分布 $q (x)$ 随机抽样得到样本 $x^*$ ，再按照均匀分布在 $(0, 1)$ 范围内抽样得到 $u$ 。
(3) 如果 $\leq \frac{p(x^*)}{cq(x^*)}$ ，则将 $x^*$ 作为抽样结果；否则，回到步骤 (2) 。
(4) 直至得到 $n$ 个随机样本，结束。

1.2 数字期望值计

一般的蒙特卡罗法，如直接抽样法、接受-拒绝抽样法、重要性抽样法，也可以用于数学期望估计 (estimation of mathematical expectation) 。假设有随机变量 $x$ ，其概率密度函数为 $p (x)$ ， $f (x)$ 为定义的函数，目标是求函数 $f (x)$ 关于密度函数 $p (x)$ 的数学期望 $E_{p(x)} [f(x)]$ 。

当 $n$ 足够大时有：
$E_{p(x)} [f(x)]\thickapprox \frac{1}{n}\sum_{i=1}^nf(x_i)$

1.3 积分计算

一般的蒙特卡罗法也可以用于定积分的近似计算，称为蒙特卡罗积分 (MonteCarlo integration) 。假设有一个函数 $h (x)$ ，如果能够将函数 $h (x)$ 分解成一个函数 $f (x)$ 和一个概率密度函数 $p (x)$ 的乘积的形式，那么就有
$\int_{\chi}h(x)dx=\int_{\chi}f(x)p(x)dx=E_{p(x)} [f(x)]\thickapprox \frac{1}{n}\sum_{i=1}^nf(x_i)$

2 马尔可夫链

2.1 基本定义

马尔可夫链：考虑一个随机变量的序列 $X=\{X_0,X_1,...,X_t,...\}$ ，这里 $X_t$ 表示时刻 $t$ 的随机变量， $t = 0, 1, 2, ...$ 。每个随机变量 $X_t(t=0,1,2,...)$ 的取值集合是相同的，称为状态空间，表示为 $S$ 。随机变量可以是离散的，也可以是连续的。以上随机变量的序列构成随机过程。

假设在时刻 $0$ 的随机变量 $X_0$ 遵循概率分布 $P(X_0)=\pi_0$ ，称为初始状态分布。在某个时刻 $t \geq 1$ 的随机变量 $X_t$ 与前一个时刻的随机变量 $X_{t−1}$ 之间有条件分布 $P(X_t|X_{t−1})$ ，如果 $X_t$ 只依赖于 $X_{t−1}$ ，而不依赖于过去的随机变量 ${X_0,X_1,...,X_{t−2}\}$ ，这一性质称为马尔可夫性，即 $P({X_t|X_0,X_1,...,X_{t−1}})=P({X_t|X_{t−1}}),t=1,2,...$ 。具有马尔可夫性的随机序列 $X=\{X_0,X_1,...,X_t,...\}$ 称为马尔可夫链，或马尔可夫过程。条件概率分布 $P({X_t|X_{t−1}})$ 称为马尔可夫链的转移概率分布。转移概率分布决定了马尔可夫链的特性。直观解释为未来只依赖于现在（假设现在已知），而与过去无关

平稳分布：设有马尔可夫链 $X=\{X_0,X_1,...,X_t,...\}$ ，其状态空间为 $S$ ，转移概率矩阵 $P=(p_{ij})$ ，如果在状态空间 $S$ 上的一个分布 $\pi = \left[ \begin{matrix} \pi_1 \\\ \pi_2 \\...\end{matrix}\right]$ 使得 $\pi=P\pi$ ，则称 $\pi$ 为马尔可夫链 $X=\{X_0,X_1,...,X_t,...\}$ 的平稳分布。

2.2 马尔可夫链的性质

马尔可夫链的性质：

①不可约。

设有马尔可夫链 $X=\{X_0,X_1,...,X_t,...\}$ ，其状态空间为 $S$ ，对于任意状态 $\in S$ ，如果存在一个时刻 $t (t > 0)$ 满足 $P(X_t=i|X_0=j)>0$ ，也就是说，时刻 $0$ 从状态 $j$ 出发，时刻 $t$ 到达状态 $i$ 的概率大于 0 ，则称此马尔可夫链 $X$ 是不可约的( irreducible ) ，否则称马尔可夫链是可约的( reducible )。

一个不可约的马尔可夫链，从任意状态出发，当经过充分长时间后，可以到达任意状态。
②非周期。

设有马尔可夫链 $X=\{X_0,X_1,...,X_t,...\}$ ，其状态空间为 $S$ ，对于任意状态 $\in S$ ，如果时刻 $0$ 从状态 $i$ 出发， $t$ 时刻返回状态 $i$ 的所有时间长 ${t:P(X_t=i|X_0=i)>0\}$ 的最大公约数是 1 ，则称此马尔可夫链 $X$ 是非周期的( aperiodic ) ，否则称马尔可夫链是周期的 (periodic) 。

一个非周期性的马尔可夫链，不存在一个状态，从这一个状态出发，再返回到这个状态时所经历的时间长呈一定的周期性。
③正常返。

设有马尔可夫链 $X=\{X_0,X_1,...,X_t,...\}$ ，其状态空间为 $S$ ，对于任意状态 $\in S$ ，定义概率 $p_{ij}^t$ 为时刻0从状态 $j$ 出发，时刻 $t$ 首次转移到状态 $i$ 的概率，即 $p_{ij}^t=P(X_t=i,X_s \neq i,s=1,2,...,t-1|X_0=j),t=1,2,...,$ ，对所有状态 $\in S$ 都满足 $\lim \limits_{t \rightarrow ∞}>0$ ，则称马尔可夫链 $X$ 是正常返的

一个正常返的马尔可夫链，其中任意一个状态，从其他任意一个状态出发，当时间趋于无穷时，首次转移到这个状态的概率不为0。

定理：不可约、非周期且正常返的马尔可夫链，有唯一平稳分布存在。
④遍历定理。

设有马尔可夫链 $X=\{X_0,X_1,...,X_t,...\}$ ，其状态空间为 $S$ ，若马尔可夫链 $X$ 是不可约、非周期且正常返的，则该马尔可夫链有唯一平稳分布 $\pi = (\pi_1,\pi_2,...,)^T$ ，并且转移概率的极限分布是马尔可夫链的平稳分布
$\lim \limits_{t \rightarrow ∞}P(X_t=i|X_0)=\pi_i,i=1,2,...,j=1,2,...$
若 $f (X)$ 是定义在状态空间上的函数， $E_{\pi}[|f(X)|]<∞$ ，则
$P(\hat{f}_t \rightarrow E_{\pi}[f(X)])=1$
，其中 $\hat{f}_t=\frac{1}{t}\sum_{s=1}^tf(x_s)$ ， $E_{\pi}[f(X)]=\sum_if(i) \pi_i$ 是 $f (X)$ 关于平稳分布 $\pi = (\pi_1,\pi_2,...,)^T$ 的的数学期望，上式表示 $\hat{f}_t \rightarrow E_{\pi}[f(X)],t \rightarrow ∞$ 几乎处处成立或以概率 1 成立。

遍历定理的直观解释：满足相应条件的马尔可夫链，当时间趋于无穷时，马尔可夫链的状态分布趋近于平稳分布，随机变量的函数的样本均值以概率1收敛于该函数的数学期望，样本均值可以认为是时间均值，而数学期望是空间均值。遍历定理实际表述了遍历性的含义：当时间趋于无穷时，时间均值等于空间均值。遍历定理的三个条件：不可约、非周期、正常返，保证了当时间趋于无穷时达到任意一个状态的概率不为 0 。

理论上并不知道经过多少次迭代，马尔可夫链的状态分布才能接近于平稳分布，在实际应用遍历定理时，取一个足够大的整数 $m$ ，经过 $m$ 次迭代之后认为状态分布就是平稳分布，这时计算从第 $m + 1$ 次迭代到第 $n$ 次迭代的均值，即 $\hat{E}f=\frac{1}{n-m}\sum_{i=m+1}^nf(x_i)$ 称为遍历均值。
⑤可逆马尔可夫链。

设有马尔可夫链 $X=\{X_0,X_1,...,X_t,...\}$ ，其状态空间为 $S$ ，庄毅概率矩阵为 $P$ ,如果有状态分布 $\pi = (\pi_1,\pi_2,...,)^T$ ，对于任意状态 $\in S$ ，对任意一个时刻 $t$ 满足 $p_{ji}\pi_j=p_{ij}\pi_i,i,j=1,2,...$ 则称此马尔可夫链 $X$ 为可逆马尔可夫链，式上式称为
细致平衡方程。

如果有可逆的马尔可夫链，那么以该马尔可夫链的平稳分布作为初始分布，进行随机状态转移，无论是面向未来还是面向过去，任何一个时刻的状态分布都是该平稳分布

3 马尔科夫链蒙特卡罗法

马尔可夫链蒙特卡罗法是以马尔可夫链为概率模型的蒙特卡罗法，它构建一个马尔可夫链，使其平稳分布就是要进行抽样的分布，首先基于该马尔可夫链进行随机游走，产生样本的序列，之后使用该平稳分布的样本进行近似数值计算。

马尔可夫链蒙特卡罗法的假设目标是对一个概率分布进行随机抽样，或者是求函数关于该概率分布的数学期望。可以采用传统的蒙特卡罗法，如接受-拒绝法、重要性抽样法，也可以使用马尔可夫链蒙特卡罗法。马尔可夫链蒙特卡罗法更适用于随机变量是多元的、密度函数是非标准形式的、随机变量各分量不独立等情况。常见的马尔可夫链蒙特卡罗法有Metropolis-Hastings算法、吉布斯抽样。

马尔可夫链蒙特卡罗法基本步骤：

①首先，在随机变量 $x$ 的状态空间 $S$ 上构造一个满足遍历定理的马尔可夫链，使其平稳分布为目标分布 $p (x)$ ；
②从状态空间的某一点 $x_0$ 出发，用构造的马尔可夫链进行随机游走，产生样本序列 $x_0,x_1,...,x_t,...$ ；
③应用马尔可夫链的遍历定理，确定正整数 $m$ 和 $n （ m < n ）$ ，得到样本集合 ${x_{m+1},x_{m+2},...,x_n\}$ ，求得函数 $f (x)$ 的均值（遍历均值） $\hat{E}f=\frac{1}{n-m}\sum_{i=m+1}^nf(x_i)$ 就是马尔可夫链蒙特卡罗法的计算公式。

4 Metropolis-Hastings 算法

Metropolis-Hastings算法是最基本的马尔可夫链蒙特卡罗法。

假设目标是对概率分布 $p (x)$ 进行抽样，构造建议分布 $q (x, x')$ ，定义接受分布 $\alpha(x,x′)$ 。进行随机游走，假设当前处于状态 $x$ ，按照建议分布 $q (x, x')$ 随机抽样，按照概率 $\alpha(x,x′)$ 接受抽样，转移到状态 $x'$ ，按照概率 $1−\alpha(x,x′)$ 拒绝抽样，停留在状态 $x$ ，持续以上操作，得到一系列样本，这样的随机游走是根据转移核为 $p(x,x′)=q(x,x′)\alpha(x,x′)$ 的可逆马尔可夫链（满足遍历定理条件）进行的，其平稳分布就是要抽样的目标分布 $p (x)$ 。

Metropolis-Hastings 算法具体如下：

输入：抽样的目标分布的密度函数 $p (x)$ ，函数 $f (x)$ ;
输出: $p (x)$ 的随机样本 $x_{m+1},x_{m+2},...,x_n$ ，函数样本均值 $f_{mn}$ .
参数：收敛步数 $m$ ，是代步数 $n$ 。
(1) 任意选择一个初始值 $x_0$
(2) 对 $i = 1, 2, \dots, n$ 循环执行
- (a) 设状态 $x_{i-1} = x$ ，按照建议分布 $q (x, x^{'})$ 随机抽取一个候选状态 $x^{'}$ 。
- (b) 计算接受概率
  $\alpha(x,x')=\min\left\{1, \frac{p(x')q(x',x)}{p(x)q(x,x')}\right \}$
- © 从区间 (0 , 1) 中按均匀分布随机抽取一个数 $u$ 。若 $\leq α (x, x')$ ，则状态 $x_i = x'$ ; 否则，状态 $x_i =x$ 。
(3) 得到样本集合 ${x_{m+1},x_{m+2},...,x_n\}$ ，计算
$f_{mn}=\frac{1}{n-m}\sum_{i=m+1}^nf(x_i)$

5 吉布斯抽样

吉布斯抽样用于多元联合分布的抽样和估计。吉布斯抽样是单分量Metropolis-Hastings算法的特殊情况。这时建议分布为满条件概率分布 $q(x,x′)=p(x'_j|x_{−j})$ ，吉布斯抽样对每次抽样的结果都接受，没有拒绝，这一点和一般的 Metropolis-Hastings 算法不同。

吉布斯抽样的基本做法是，从联合分布定义满条件概率分布，依次从满条件概率分布进行抽样，得到联合分布的随机样本。假设多元联合概率分布为 $p(x)=p(x_1,x_2,...,x_k)$ ，吉布斯抽样从一个初始样本 $x^{(0)}=(x_1^{(0)},x_2^{(0)},...,x_k^{(0)})^T$ 出发，不断进行迭代，每一次迭代得到联合分布的一个样本 $x^{(i)}=(x_1^{(i)},x_2^{(i)},...,x_k^{(i)})^T$ 。在第 $i$ 次迭代中，依次对第 $j$ 个变量按照满条件概率分布随机抽样， $p(x_j|x_1^{(i)},x_2^{(i)},...,x_{j-1}{(i)},x_{j+1}{(i-1)} ,x_{k}{(i-1)}) ,j=1,2,...,k$ ，得到 $x_j^{(i)}$ 。最终得到样本序列 ${x^{(0)},x^{(1)},...,x^{(n)}\}$

吉布斯抽样算法具体如下：

输入：抽样的目标分布的密度函数 $p (x)$ ，函数 $f (x)$ ;
输出: $p (x)$ 的随机样本 $x_{m+1},x_{m+2},...,x_n$ ，函数样本均值 $f_{mn}$ .
参数：收敛步数 $m$ ，是代步数 $n$ 。
（1）初始化。给出初始样本 $x^{(0)}=(x_1^{(0)},x_2^{(0)},...,x_k^{(0)})^T$
（2）对 $i$ 进行循环执行。设第 $(i - 1)$ 次迭代结束的样本为 $x^{(i-1)}=(x_1^{(i-1)},x_2^{(i-1)},...,x_k^{(i-1)})^T$ ，则第 $i$ 次迭代进行如下几步操作
- (1) 由满条件分布 $p(x_1|x_2^{(i-1)},... ,x_{k}^{(i-1)})$ 抽取 $x_1^{(i)}$
- ( $j$ ) 由满条件分布 $p(x_j|x_1^{(i)},x_2^{(i)},...,x_{j-1}^{(i)},x_{j+1}^{(i-1)} ,x_{k}^{(i-1)})$ 抽取 $x_j^{(i)}$
- ( $k$ ) 由满条件分布 $p(x_k|x_1^{(i)},... ,x_{k-1}^{(i)})$ 抽取 $x_k^{(i)}$
得到第 $i$ 次迭代值 $x^{(i)}=(x_1^{(i)},x_2^{(i)},...,x_k^{(i)})^T$
(3) 得到样本集合 ${x_{m+1},x_{m+2},...,x_n\}$
(4) 计算
$f_{mn}=\frac{1}{n-m}\sum_{i=m+1}^nf(x^{(i)})$

healed萌

关注

23
点赞
踩
18

收藏

觉得还不错? 一键收藏
打赏
0
评论
29、马尔可夫链蒙特卡罗法

本文介绍了马尔可夫链蒙特卡罗法，首先介绍了蒙特卡罗法、马尔可夫链的概念，然后进行介绍马尔可夫链蒙特卡罗法，最后介绍了两种常用的马尔可夫链蒙特卡罗法（Metropolis-Hastings算法、吉布斯抽样）
复制链接

扫一扫