LDA-math-MCMC 和 Gibbs Sampling
来源:互联网 发布:淘宝的聚划算怎么用 编辑:程序博客网 时间:2024/04/29 22:23
随机模拟
随机模拟(或者统计模拟)方法有一个很酷的别名是蒙特卡罗方法(Monte Carlo Simulation)。这个方法的发展始于20世纪40年代,和原子弹制造的曼哈顿计划密切相关,当时的几个大牛,包括乌拉姆、冯.诺依曼、费米、费曼、Nicholas Metropolis, 在美国洛斯阿拉莫斯国家实验室研究裂变物质的中子连锁反应的时候,开始使用统计模拟的方法,并在最早的计算机上进行编程实现。
随机模拟与计算机
现代的统计模拟方法最早由数学家乌拉姆提出,被Metropolis命名为蒙特卡罗方法,蒙特卡罗是著名的赌场,赌博总是和统计密切关联的,所以这个命名风趣而贴切,很快被大家广泛接受。被不过据说费米之前就已经在实验中使用了,但是没有发表。说起蒙特卡罗方法的源头,可以追溯到18世纪,布丰当年用于计算
蒙特卡罗方法
统计模拟中有一个重要的问题就是给定一个概率分布
生成一个概率分布的样本
而我们常见的概率分布,无论是连续的还是离散的分布,都可以基于
[Box-Muller 变换] 如果随机变量
则
其它几个著名的连续分布,包括指数分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等,也都可以通过类似的数学变换得到;离散的分布通过均匀分布更加容易生成。更多的统计分布如何通过均匀分布的变换生成出来,大家可以参考统计计算的书,其中 Sheldon M. Ross 的《统计模拟》是写得非常通俗易懂的一本。
不过我们并不是总是这么幸运的,当
p(x)=p˜(x)∫p˜(x)dx ,而p˜(x) 我们是可以计算的,但是底下的积分式无法显式计算。p(x,y) 是一个二维的分布函数,这个函数本身计算很困难,但是条件分布p(x|y),p(y|x) 的计算相对简单;如果p(x) 是高维的,这种情形就更加明显。
此时就需要使用一些更加复杂的随机模拟的方法来生成样本。而本节中将要重点介绍的 MCMC(Markov Chain Monte Carlo) 和 Gibbs Sampling算法就是最常用的一种,这两个方法在现代贝叶斯分析中被广泛使用。要了解这两个算法,我们首先要对马氏链的平稳分布的性质有基本的认识。
3.2 马氏链及其平稳分布
马氏链的数学定义很简单
也就是状态转移的概率只依赖于前一个状态。
我们先来看马氏链的一个具体的例子。社会学家经常把人按其经济状况分成3类:下层(lower-class)、中层(middle-class)、上层(upper-class),我们用1,2,3 分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别,那么他的孩子属于下层收入的概率是 0.65, 属于中层收入的概率是 0.28, 属于上层收入的概率是 0.07。事实上,从父代到子代,收入阶层的变化的转移概率如下
使用矩阵的表示方式,转移概率矩阵记为
假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量
假设初始概率分布为
我们发现从第7代人开始,这个分布就稳定不变了,这个是偶然的吗?我们换一个初始概率分布
我们发现,到第9代人的时候, 分布又收敛了。最为奇特的是,两次给定不同的初始概率分布,最终都收敛到概率分布
我们发现,当
马氏链定理: 如果一个非周期马氏链具有转移概率矩阵
limn→∞Pn=⎡⎣⎢⎢⎢⎢⎢π(1)π(1)⋯π(1)⋯π(2)π(2)⋯π(2)⋯⋯⋯⋯⋯⋯π(j)π(j)⋯π(j)⋯⋯⋯⋯⋯⋯⎤⎦⎥⎥⎥⎥⎥
-
π(j)=∑i=0∞π(i)Pij -
π 是方程πP=π 的唯一非负解
其中,
这个马氏链的收敛定理非常重要,所有的 MCMC(Markov Chain Monte Carlo) 方法都是以这个定理作为理论基础的。 定理的证明相对复杂,一般的随机过程课本中也不给证明,所以我们就不用纠结它的证明了,直接用这个定理的结论就好了。我们对这个定理的内容做一些解释说明:
- 该定理中马氏链的状态不要求有限,可以是有无穷多个的;
- 定理中的“非周期“这个概念我们不打算解释了,因为我们遇到的绝大多数马氏链都是非周期的;
- 两个状态
i,j 是连通并非指i 可以直接一步转移到j (Pij>0 ),而是指i 可以通过有限的n 步转移到达j (Pnij>0 )。马氏链的任何两个状态是连通的含义是指存在一个n , 使得矩阵Pn 中的任何一个元素的数值都大于零。 - 我们用
Xi 表示在马氏链上跳转第i 步后所处的状态,如果limn→∞Pnij=π(j) 存在,很容易证明以上定理的第二个结论。由于P(Xn+1=j)=∑i=0∞P(Xn=i)P(Xn+1=j|Xn=i)=∑i=0∞P(Xn=i)Pij
上式两边取极限就得到π(j)=∑i=0∞π(i)Pij
从初始概率分布
由马氏链收敛的定理, 概率分布
所以
3.2 Markov Chain Monte Carlo
对于给定的概率分布
这个绝妙的想法在1953年被 Metropolis想到了,为了研究粒子系统的平稳性质, Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题,首次提出了基于马氏链的蒙特卡罗方法,即Metropolis算法,并在最早的计算机上编程实现。Metropolis 算法是首个普适的采样方法,并启发了一系列 MCMC方法,所以人们把它视为随机模拟技术腾飞的起点。 Metropolis的这篇论文被收录在《统计学中的重大突破》中, Metropolis算法也被遴选为二十世纪的十个最重要的算法之一。
我们接下来介绍的MCMC 算法是 Metropolis 算法的一个改进变种,即常用的 Metropolis-Hastings 算法。由上一节的例子和定理我们看到了,马氏链的收敛性质主要由转移矩阵
定理:[细致平稳条件] 如果非周期马氏链的转移矩阵
则
其实这个定理是显而易见的,因为细致平稳条件的物理含义就是对于任何两个状态
由于
假设我们已经有一个转移矩阵为
取什么样的
于是(*)式就成立了。所以有
于是我们把原来具有转移矩阵
在改造
马氏链转移和接受概率
假设我们已经有一个转移矩阵Q(对应元素为
上述过程中
以上的 MCMC 采样算法已经能很漂亮的工作了,不过它有一个小的问题:马氏链
假设
上式两边扩大5倍,我们改写为
看,我们提高了接受率,而细致平稳条件并没有打破!这启发我们可以把细致平稳条件(**) 式中的
于是,经过对上述MCMC 采样算法中接受率的微小改造,我们就得到了如下教科书中最常见的 Metropolis-Hastings 算法。
对于分布
此处
那么以上的 Metropolis-Hastings 算法一样有效。
3.2 Gibbs Sampling
对于高维的情形,由于接受率
所以得到
即
基于以上等式,我们发现,在
平面上马氏链转移矩阵的构造
于是我们可以如下构造平面上任意两点之间的转移概率矩阵Q
有了如上的转移矩阵 Q, 我们很容易验证对平面上任意两点
于是这个二维空间上的马氏链将收敛到平稳分布
二维Gibbs Sampling 算法中的马氏链转移
以上采样过程中,如图所示,马氏链的转移只是轮换的沿着坐标轴
以上的过程我们很容易推广到高维的情形,对于(***) 式,如果
此时转移矩阵 Q 由条件分布
- 如果当前状态为
(x1,x2,⋯,xn) ,马氏链转移的过程中,只能沿着坐标轴做转移。沿着xi 这根坐标轴做转移的时候,转移概率由条件概率p(xi|x1,⋯,xi−1,xi+1,⋯,xn) 定义; - 其它无法沿着单根坐标轴进行的跳转,转移概率都设置为 0。
于是我们可以把Gibbs Smapling 算法从采样二维的
以上算法收敛后,得到的就是概率分布
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling
- LDA-math-MCMC 和 Gibbs Sampling (我爱NLP)
- 马尔科夫链MCMC采样算法和LDA Gibbs Sampling
- LDA文本建模(2)——MCMC和Gibbs Sampling
- MCMC和ÚGibbs Sampling
- MCMC与Gibbs Sampling
- Gibbs Sampling实现LDA
- Lda gibbs sampling --- python
- applicationContext.getBean() 随处可获取basedao
- BeegoAdmin 使用总结
- magento mysql 配置
- 金蝶委外入库核销慢
- 【ACM之旅】数列排序
- LDA-math-MCMC 和 Gibbs Sampling
- BMFont中文字体图集制作的方法~(for unity ngui)
- Python个人代码库 rom里提取apk脚本
- Android设计模式源码解析之ListView观察者模式
- 最好的5个Android ORM框架
- JVM最多能创建多少个线程: unable to create new native thread
- RTP协议分析
- Real-time Compressive Tracking
- linux IPC之信号