随机模拟中有一个重要的问题就是给定一个概率分布p(x), 如何在计算机中生成其样本。
常见的分布,如正态分布、指数分布、Gamma分布、t分布等都可以通过均匀分布Uniform(0, 1)变换得到,而均匀分布可通过线性同余发生器生成。例如,正态分布可由Box-Muller变化而来:
其中U1和U2相互独立且服从均匀分布,则Z0和Z1也相互独立且服从标准正态分布。
但当概率分布p(x)形式很复杂或是个高维分布时,样本的生成很困难。
在此,我们介绍两种随机模拟的方法来生成样本:MCMC(Markov Chain Monte Carlo) 和 Gibbs Sampling。
要了解这两个算法,首先要对马尔可夫链(Markov Chain)以及马尔可夫平稳分布要有基本认识。
马尔可夫链的定义为:
即随机向量(X)转移的概率只依赖于前一个随机向量。
一个具体的例子:
社会学家经常把人按其经济状况分成 3 类:下层 (lower-class)、中层 (middle-class)、上层 (upper-class),我们用 1、2、3 分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别,那么他的孩子属于下层收入的概率是 0.65,属于中层收入的概率是 0.28,属于上层收入的概率是 0.07。事实上,从父代到子代,收入阶层的变化的转移概率如下:
使用矩阵的表示方式,转移概率矩阵记为:
假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量 ,那么他们的子女的分布比例将是
他们的孙子代的分布比例将是
第n代子孙的收入分布比例是
给定不同的初始概率分布,我们发现从第七代开始都收敛到
也就是说,最后收敛的概率分布与初始概率分布无关,只与转移矩阵相关。
马尔可夫链定理: 如果一个非周期马尔可夫链具有转移概率矩阵$P$, 且它的任何两个状态是连通的,那么 存在且与i无关
π称为马尔可夫链的平稳分布。
MCMC
对于给定的概率分布p(x),如果我们能构造一个转移矩阵$P$的马尔可夫链,使其平稳分布正好是p(x)。
所以基于马尔可夫链做采样的关键问题是如何构造转移矩阵$P$,我们主要需要使用细致平稳条件定理:
其证明过程如下:
但如果马尔可夫链还未达到平稳状态,此时:
其中q(i, j)表示从状态i转移到状态j的概率。
因此我们引入a(i, j) = p(j)q(j, i); a(j, i) = p(i)q(j, i),使得:
此时,具有Q'的马尔可夫链的平稳分布就是p(x)。