随机模拟（MCMC）

首页 > 代码库 > 随机模拟（MCMC）

2024-10-13 13:09:39 214人阅读

http://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/

http://blog.csdn.net/lin360580306/article/details/51240398

随机模拟(或者统计模拟)方法有一个很酷的别名是蒙特卡罗方法(Monte Carlo Simulation)。这个方法的发展始于20世纪40年代，和原子弹制造的曼哈顿计划密切相关，当时的几个大牛，包括乌拉姆、冯.诺依曼、费米、费曼、Nicholas Metropolis，在美国洛斯阿拉莫斯国家实验室研究裂变物质的中子连锁反应的时候，开始使用统计模拟的方法,并在最早的计算机上进行编程实现。

随机模拟与计算机

现代的统计模拟方法最早由数学家乌拉姆提出，被Metropolis命名为蒙特卡罗方法，蒙特卡罗是著名的赌场，赌博总是和统计密切关联的，所以这个命名风趣而贴切，很快被大家广泛接受。被不过据说费米之前就已经在实验中使用了，但是没有发表。说起蒙特卡罗方法的源头，可以追溯到18世纪，布丰当年用于计算 $π$

蒙特卡罗方法

统计模拟中有一个重要的问题就是给定一个概率分布 $p (x)$

技术分享

生成一个概率分布的样本

而我们常见的概率分布，无论是连续的还是离散的分布，都可以基于 $U n i f o r m (0, 1)$

[Box-Muller 变换] 如果随机变量 $U_{1}, U_{2}$

\begin{aligned} Z_{0} & = \sqrt{- 2 \ln U_{1}} c o s (2 π U_{2}) \\ Z_{1} & = \sqrt{- 2 \ln U_{1}} s i n (2 π U_{2}) \end{aligned}

则 $Z_{0}, Z_{1}$

其它几个著名的连续分布，包括指数分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等,也都可以通过类似的数学变换得到；离散的分布通过均匀分布更加容易生成。更多的统计分布如何通过均匀分布的变换生成出来，大家可以参考统计计算的书，其中 Sheldon M. Ross 的《统计模拟》是写得非常通俗易懂的一本。

不过我们并不是总是这么幸运的，当 $p (x)$

$p (x) = \frac{\tilde{p} (x)}{\int \tilde{p} (x) d x}$
$p (x, y)$

此时就需要使用一些更加复杂的随机模拟的方法来生成样本。而本节中将要重点介绍的 MCMC(Markov Chain Monte Carlo) 和 Gibbs Sampling算法就是最常用的一种，这两个方法在现代贝叶斯分析中被广泛使用。要了解这两个算法，我们首先要对马氏链的平稳分布的性质有基本的认识。

3.2 马氏链及其平稳分布

马氏链的数学定义很简单

P (X_{t + 1} = x | X_{t}, X_{t - 1}, \dots) = P (X_{t + 1} = x | X_{t})

$Z_{0}, Z_{1}$

我们先来看马氏链的一个具体的例子。社会学家经常把人按其经济状况分成3类：下层(lower-class)、中层(middle-class)、上层(upper-class)，我们用1,2,3 分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是 0.65, 属于中层收入的概率是 0.28, 属于上层收入的概率是 0.07。事实上，从父代到子代，收入阶层的变化的转移概率如下

技术分享

使用矩阵的表示方式，转移概率矩阵记为

P = [\begin{matrix} 0.65 & 0.28 & 0.07 \\ 0.15 & 0.67 & 0.18 \\ 0.12 & 0.36 & 0.52 \end{matrix}]

假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量 $π_{0} = [π_{0} (1), π_{0} (2), π_{0} (3)]$

假设初始概率分布为 $π_{0} = [0.21, 0.68, 0.11]$

技术分享

我们发现从第7代人开始，这个分布就稳定不变了，这个是偶然的吗？我们换一个初始概率分布 $π_{0} = [0.75, 0.15, 0.1]$

技术分享

我们发现，到第9代人的时候, 分布又收敛了。最为奇特的是，两次给定不同的初始概率分布，最终都收敛到概率分布 $π = [0.286, 0.489, 0.225]$

P^{20} = P^{21} = \dots = P^{100} = \dots = [\begin{matrix} 0.286 & 0.489 & 0.225 \\ 0.286 & 0.489 & 0.225 \\ 0.286 & 0.489 & 0.225 \end{matrix}]

我们发现，当 $n$

马氏链定理： 如果一个非周期马氏链具有转移概率矩阵 $P$

$lim_{n \to \infty} P^{n} = [\begin{matrix} π (1) & π (2) & \dots & π (j) & \dots \\ π (1) & π (2) & \dots & π (j) & \dots \\ \dots & \dots & \dots & \dots & \dots \\ π (1) & π (2) & \dots & π (j) & \dots \\ \dots & \dots & \dots & \dots & \dots \end{matrix}]$

$π (j) = \sum_{i = 0}^{\infty} π (i) P_{i j}$
$π$

其中,

π = [π (1), π (2), \dots, π (j), \dots], \sum_{i = 0}^{\infty} π_{i} = 1

$Z_{0}, Z_{1}$

这个马氏链的收敛定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以这个定理作为理论基础的。 定理的证明相对复杂，一般的随机过程课本中也不给证明，所以我们就不用纠结它的证明了，直接用这个定理的结论就好了。我们对这个定理的内容做一些解释说明：

该定理中马氏链的状态不要求有限，可以是有无穷多个的；
定理中的“非周期“这个概念我们不打算解释了，因为我们遇到的绝大多数马氏链都是非周期的；
两个状态 $i, j$
我们用 $X_{i}$
$\begin{aligned} P (X_{n + 1} = j) & = \sum_{i = 0}^{\infty} P (X_{n} = i) P (X_{n + 1} = j | X_{n} = i) \\ = \sum_{i = 0}^{\infty} P (X_{n} = i) P_{i j} \end{aligned}$
$X_{i}$

从初始概率分布 $π_{0}$

\begin{aligned} X_{0} & \sim π_{0} (x) \\ X_{i} & \sim π_{i} (x), π_{i} (x) = π_{i - 1} (x) P = π_{0} (x) P^{n} \end{aligned}

$Z_{0}, Z_{1}$

\begin{aligned} X_{0} & \sim π_{0} (x) \\ X_{1} & \sim π_{1} (x) \\ \dots \\ X_{n} & \sim π_{n} (x) = π (x) \\ X_{n + 1} & \sim π (x) \\ X_{n + 2} & \sim π (x) \\ \dots \end{aligned}

$Z_{0}, Z_{1}$

3.3 Markov Chain Monte Carlo

对于给定的概率分布 $p (x)$

这个绝妙的想法在1953年被 Metropolis想到了，为了研究粒子系统的平稳性质， Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。 Metropolis的这篇论文被收录在《统计学中的重大突破》中， Metropolis算法也被遴选为二十世纪的十个最重要的算法之一。

我们接下来介绍的MCMC 算法是 Metropolis 算法的一个改进变种，即常用的 Metropolis-Hastings 算法。由上一节的例子和定理我们看到了，马氏链的收敛性质主要由转移矩阵 $P$

定理：[细致平稳条件] 如果非周期马氏链的转移矩阵 $P$

\begin{matrix} (1) & π (i) P_{i j} = π (j) P_{j i} for all i, j \end{matrix}

$Z_{0}, Z_{1}$

其实这个定理是显而易见的，因为细致平稳条件的物理含义就是对于任何两个状态 $i, j$

\begin{aligned} \sum_{i = 1}^{\infty} π (i) P_{i j} = \sum_{i = 1}^{\infty} π (j) P_{j i} = π (j) \sum_{i = 1}^{\infty} P_{j i} = π (j) \\ \Rightarrow π P = π \end{aligned}

$Z_{0}, Z_{1}$

假设我们已经有一个转移矩阵为 $Q$

p (i) q (i, j) \neq p (j) q (j, i)

$Z_{0}, Z_{1}$

\begin{matrix} (2) & p (i) q (i, j) α (i, j) = p (j) q (j, i) α (j, i) (*) \end{matrix}

$Z_{0}, Z_{1}$

α (i, j) = p (j) q (j, i) ， α (j, i) = p (i) q (i, j)

$Z_{0}, Z_{1}$

\begin{matrix} (3) & p (i) \underset{Q^{'} (i, j)}{\underset{⏟}{q (i, j) α (i, j)}} = p (j) \underset{Q^{'} (j, i)}{\underset{⏟}{q (j, i) α (j, i)}} (* *) \end{matrix}

$Z_{0}, Z_{1}$

在改造 $Q$

马氏链转移和接受概率

假设我们已经有一个转移矩阵Q(对应元素为 $q (i, j)$

技术分享

上述过程中 $p (x), q (x | y)$

以上的 MCMC 采样算法已经能很漂亮的工作了，不过它有一个小的问题：马氏链 $Q$

假设 $α (i, j) = 0.1, α (j, i) = 0.2$

p (i) q (i, j) \times 0.1 = p (j) q (j, i) \times 0.2

$Z_{0}, Z_{1}$

p (i) q (i, j) \times 0.5 = p (j) q (j, i) \times 1

$Z_{0}, Z_{1}$

α (i, j) = min {\frac{p (j) q (j, i)}{p (i) q (i, j)}, 1}

$Z_{0}, Z_{1}$

技术分享

对于分布 $p (x)$

p (x) Q^{'} (x \to y) = p (y) Q^{'} (y \to x)

$Z_{0}, Z_{1}$

p (x) Q^{'} (x \to y) = p (y) Q^{'} (y \to x)

$Z_{0}, Z_{1}$

3.2 Gibbs Sampling

对于高维的情形，由于接受率 $α$

\begin{aligned} p (x_{1}, y_{1}) p (y_{2} | x_{1}) = p (x_{1}) p (y_{1} | x_{1}) p (y_{2} | x_{1}) \\ p (x_{1}, y_{2}) p (y_{1} | x_{1}) = p (x_{1}) p (y_{2} | x_{1}) p (y_{1} | x_{1}) \end{aligned}

$Z_{0}, Z_{1}$

\begin{matrix} (4) & p (x_{1}, y_{1}) p (y_{2} | x_{1}) = p (x_{1}, y_{2}) p (y_{1} | x_{1}) (* * *) \end{matrix}

$Z_{0}, Z_{1}$

p (A) p (y_{2} | x_{1}) = p (B) p (y_{1} | x_{1})

$Z_{0}, Z_{1}$

p (A) p (x_{2} | y_{1}) = p (C) p (x_{1} | y_{1}) .

技术分享 平面上马氏链转移矩阵的构造

于是我们可以如下构造平面上任意两点之间的转移概率矩阵Q

\begin{aligned} Q (A \to B) & = p (y_{B} | x_{1}) & 如果 x_{A} = x_{B} = x_{1} \\ Q (A \to C) & = p (x_{C} | y_{1}) & 如果 y_{A} = y_{C} = y_{1} \\ Q (A \to D) & = 0 & 其它 \end{aligned}

有了如上的转移矩阵 Q, 我们很容易验证对平面上任意两点 $X, Y$

p (X) Q (X \to Y) = p (Y) Q (Y \to X)

$Z_{0}, Z_{1}$

技术分享

Gibbs Sampling 算法中的马氏链转移

以上采样过程中，如图所示，马氏链的转移只是轮换的沿着坐标轴 $x$

以上的过程我们很容易推广到高维的情形，对于(***) 式，如果 $x_{1}$

\begin{matrix} (5) & p (x_{1}, y_{1}) p (y_{2} | x_{1}) = p (x_{1}, y_{2}) p (y_{1} | x_{1}) \end{matrix}

$Z_{0}, Z_{1}$

如果当前状态为 $(x_{1}, x_{2}, \dots, x_{n})$
其它无法沿着单根坐标轴进行的跳转，转移概率都设置为 0。

于是我们可以把Gibbs Smapling 算法从采样二维的 $p (x, y)$

技术分享

以上算法收敛后，得到的就是概率分布 $p (x_{1}, x_{2}, \dots, x_{n})$

随机模拟（MCMC）

声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉：投诉/举报工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。

联系
我们

首页 > 代码库 > 随机模拟（MCMC）

随机模拟（MCMC）

看完仍有疑问？有类似问题直接问程序猿