Gibbs Sampling(一)：随机数产生方法介绍 & Monte Carlo方法

转自：http://blog.sina.com.cn/s/blog_5033f3b40101jfqu.html

本系列是针对gibbs采样在LDA中的应用而学习总结的四个部分，但并不着重强调其在LDA中的应用，而是作为一个应用的普通例子来general的介绍Gibbs sampling的前前后后，主要参考是Sheldon M.Ross的《Simulation》（Fourth Edition）边学习边补充边总结，有不足或者理解不对的敬请指出.

一、如何产生伪随机数？或者说，如何产生0,1之间均匀分布的随机数？(1.产生方法；2.Monte Carlo方法)(注：这里是后续过程的基础，基本上所有的服从某一分布的随机数的产生方法都可以或者会用到uniform的random variables)

二、如何产生服从某一离散分布的样本点？(1.general的方法；2.example1-计算均值；3.如何产生伯努利分布的随机变量；4.如何产生泊松分布的随机变量；5.产生二项分布的随机变量)(注：这里因为牵涉到在lda的gibbs采样实现过程中对多项式分布采样，所以需要了解这一点)

三、如何产生服从某一连续分布的样本点？

四、MCMC & Gibbs Sampling 介绍

暂定为四个部分，其中，一为基础知识，二三是采样中经常会用到的方法，在通常的文章中都是一句话带过，但是每次都不懂细节怎么实现，因此也一并在此总结介绍，但就算不理解他们并不会影响对Gibbs采样的大致理解。

==========================================================

后补：现在已经基本完成，但是分了三个部分：如下

关于Gibbs Sampling & LDA笔记一：随机数产生方法介绍 & Monte Carlo方法

关于Gibbs Sampling & LDA笔记二：Gibbs Sampling总结

关于Gibbs Sampling & LDA笔记三：补充

===================================================

====================================================

一、如何产生伪随机数？(Pseudorandom Number Generation)

1.产生方法

    一种常用的产生伪随机数的方法是首先设定一个初始值x₀,我们称其为种子(seed)，然后根据下式递归的计算x_n(n>=1),

x_n = ax_n-1 modulo m

其中a和m是固定的整数，modulo是取余的意思。这个式子表明每个xn就是0,1,2…m-1中的一个，xn/m就称为一个伪随机数，并被用来作为对(0,1)之间均匀分布的随机变量的一个估计样本点。

    因为每个数字xn都一定会取0,1,…,m-1之间的值，所以当一定次数之后，所产生的x一定会有重复(最多m次就会有重复), 根据上面的递归公式，当某点开始重复之后，那么这个点之后的样本点也都是对之前的重复。因此我们的目的就是让最开始重复的这个点尽量的延后。

        一个常用的guideline就是m尽可能的大。例如，对于32-bit的machine来说，选择m=2³¹-1, a=7⁵=16807时可以取得最佳效果。

        还有其他的方法来生成随机数，但这里我们仅仅是举个例子来说明如何产生随机数。对以下的讨论，我们都是假设我们有一个black box黑盒子，这个黑盒子可以给出满足要求的随机数(0,1之间均匀分布的随机数)。

2. Monte Carlo方法

说到monte carlo方法就必须得说起积分问题，因为monte Carlo方法最初就是被设计来求解积分的，这也是随机数（random numbers）的最早应用。

【插入一句，积分在概率分布中很重要的应用就是求相对于后验分布的期望，比如一个模型的观测是x,我们想要求某个变量y的概率p(y|x),假设模型参数为theta,就有p(y|x)=∫p(y|θ,x)p(θ|x)dθ,很多时候p(y|θ,x)是很容易知道的，比如可以是高斯分布等，而往往p(θ|x)这个后验比较难求，或者说比较不容易写成某一个已有的单一分布的形式，这时候就可以用采样的方法，也就是从这个分布p(θ|x)中采样，(因为分布形式已知，就可以知道θ取各个值的概率，就可以采样)，如得到N个θi(i=1~N)，然后用平均值∑(p(y|θi,x))/N来代替原来的积分。】

2.1 [0,1]区间的积分

假设我们需要求解积分：

( 1 ） (编辑不了公式，只能从word截图过来的，将就将就。。)

为计算\theta的值，我们注意到，如果x是(0,1)之间的均匀分布，x ~ uniform(0,1), 我们可以把上面的积分式改写成：

（ 2 ）

通过上式，我们把这个积分转换成了期望的形式，这是引出monte carlo的重要的一步。咱们接着往下看。

    假入我们有k个独立的服从0,1之间均匀分布的随机变量U1，U2，…,Uk,那么就有g(U1), g(U2),…, g(Uk)也是独立的且同分布并且均值都是\theta. （注：这里如何得到他们就要涉及到上面第一部分讲的如何产生随机数了，这里我们假设一个完美的黑盒子来产生他们，而不去理会产生的细节过程）。根据大数定律，当k足够大时，就有下面的式子成立：

（ 3 ）

再联合2,3两个式子就发现，我们就得到

（ 4 ）

由最终的（4）式，我们就发现一个重要的事情，就是对于原来的积分，我们可以通过产生大量的随机数U_i，然后取均值来作为的估计值，这就是Monte Carlo方法(the Monte Carlo approach)。这里并不是严格的定义，just an introduction, 只是说这种应用产生的一系列随机数来进行估计的方法是Monte Carlo方法。        

2.2 任意有限区间的积分

上面我们得到的是0,1之间的积分，如果变换了积分范围，怎么办？万变不离其宗，不管怎样就是要通过变量替换将积分范围变换到0,1之间。

如对于积分, 假设新变换后处于0,1之间的变量为y，则容易的y与x的变换为y=(x-a)/(b-a), 并有dy=dx/(b-a), 代入积分就有：

这里就可以通过生成0,1之间的随机数y，然后代入h(y),求平均值就可以求出。

2.3 包含有无限积分的形式：

        假如我们要求的是,原理同上，依然是假设变换后处于[0,1]区间内的变量为y，则容易得到y与x之间的映射关系可以建立为y=1/(x+1), 从而有dy=-dx/(x+1)²=-y²dx,代入原积分式，有

总结

从2.1,2.2,2.3总结，其实用monte carlo方法(Monte Carlo approach)的关键是要能够吧积分表达成独立的随机变量的期望的形式，即：

其中U1,U2…,Un是[0,1]之间独立的均匀分布的变量。因此，我们进行k次独立实验

                                      

然后是独立同分布的且期望为,就可以用平均值来估计。

重要！！！在大部分的从某一概率分布产生样本点的过程中都会用到概率函数的累计分布，，而累积分布cdf是由pdf进行区间积分求得的，(因为累积分布是在[0,1]之间取值的，就转化到上面的问题上去了)，因此这里的各种积分情况的介绍对于后续的理解是有极大的帮助的。

3. Monte Carlo外传

ross的书里是没有详细介绍这个的，下面主要是我参考其他资料总结的MonteCarlo方法,与大家分享。

【关于Monte Carlo的八卦】 说到Monte Carlo的渊源，很多人都会意外的。这个方法"Monte Carlo"最初是为了建造更好的原子弹而设计的！！！ 吓到了吧，木有想到吧，再细想一下就可以理解了，Monte Carlo牵涉到什么，积分啊，积分可是会经常在那些高端技术里用到的，所以这就好理解了。而最初使用这个方法的那两个哥们不知道咋想的，用一个赌场的名字来作为这个方法的名字，你没看错，Monte Carlo在当初是Ulam(最初使用的两个人之一)的叔叔经常去赌博的一个赌场的名字。

【关于MonteCarlo的定义】目前为止，并没有明确地Monte Carlo的定义，有些人会用随机模拟(stochastic simulation)来指代任何可以指代的东西，而仅仅对于Monte Carlo积分和MonteCarlo测试才会使用Monte Carlo这个词，而还有些人不太care到底应该在哪些场合需要用这个词。但是我看到一个定义觉得写得挺好的："Monte Carlo is the art of approximating an expectation by the sample mean of a function of simulated random variables"，就是Monte Carlo可以认为是一种通过模拟的随机变量的函数的均值来估计某个期望的方法。我们再联系上面讲的，就会很好理解这个话了。

【引用wiki】: Monte Carlo methods (or Monte Carlo experiments) are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results; They are often used in physical and mathematical problems and are most suited to be applied when it is impossible to obtain a closed-form expression or infeasible to apply a deterministic algorithm. Monte Carlo methods are mainly used in three distinct problems: optimization, numerical integration and generation of samples from a probability distribution.

【可以作为总结的一段话：】Monte Carlo sampling is often used in two kinds of related problems.

• Sampling from a distributionp(x), often a posterior distribution.

• Computing approximate integrals of the form i.e., computing expectation of f(x) using density p(x).

The above problems are related because if we can sample from p(x) then, we can also solve the problem of computing integrals.

所以下一篇是关于如何从一个已知的分布进行采样的笔记。

参考资料：

1. Sheldon M.Ross的《Simulation》,(Fourth Edition)

2. http://omlc.ogi.edu/news/sep98/montecarlosampling/montecarlo1.html

3. Monte Carlo Methods and Importance Sampling

4. Wikipedia_'Monte Carlo methods'

欢迎各位讨论拍砖。

=====================================

本文地址：http://blog.sina.com.cn/s/blog_5033f3b40101jfqu.html

================================================