1 随机模拟

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

随机模拟与计算机

现代的统计模拟方法最早由数学家乌拉姆提出，被Metropolis命名为蒙特卡罗方法，蒙特卡罗是著名的赌场，赌博总是和统计密切关联的，所以这个命名风趣而贴切，很快被大家广泛接受。被不过据说费米之前就已经在实验中使用了，但是没有发表。说起蒙特卡罗方法的源头，可以追溯到18世纪，布丰当年用于计算 π 的著名的投针实验就是蒙特卡罗模拟实验。统计采样的方法其实数学家们很早就知道，但是在计算机出现以前，随机数生成的成本很高，所以该方法也没有实用价值。随着计算机技术在二十世纪后半叶的迅猛发展，随机模拟技术很快进入实用阶段。对那些用确定算法不可行或不可能解决的问题，蒙特卡罗方法常常为人们带来希望。

蒙特卡罗方法

统计模拟中有一个重要的问题就是给定一个概率分布 p(x) ，我们如何在计算机中生成它的样本。一般而言均匀分布  Uniform(0,1) 的样本是相对容易生成的。通过线性同余发生器可以生成伪随机数，我们用确定性算法生成 [0,1]   之间的伪随机数序列后，这些序列的各种统计指标和均匀分布  Uniform(0,1)  的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。

生成一个概率分布的样本

而我们常见的概率分布，无论是连续的还是离散的分布，都可以基于 Uniform(0,1) 的样本生成。例如正态分布可以通过著名的Box-Muller 变换得到

[Box-Muller 变换] 如果随机变量  U1,U2  独立且 U1,U2∼Uniform[0,1]，

则  Z0,Z1  独立且服从标准正态分布。

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

不过我们并不是总是这么幸运的，当 p(x) 的形式很复杂，或者  p(x)  是个高维的分布的时候，样本的生成就可能很困难了。譬如有如下的情况

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

2 马氏链及其平稳分布

马氏链的数学定义很简单

也就是状态转移的概率只依赖于前一个状态。

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

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

我们发现从第7代人开始，这个分布就稳定不变了，这个是偶然的吗？我们换一个初始概率分布 π0 =[0.75,0.15,0.1] 试试看，继续计算前 n 代人的分布状况如下

我们发现，到第9代人的时候,分布又收敛了。最为奇特的是，两次给定不同的初始概率分布，最终都收敛到概率分布  π=[0.286,0.489,0.225] ，也就是说收敛的行为和初始概率分布  π0  无关。这说明这个收敛行为主要是由概率转移矩阵 P 决定的。

（接二）