Gibbs sampling

　Gibbs sampling(吉布斯采样)(资料集合)

　　维基百科，自由的百科全书:

　　In statistics, Gibbs sampling or a Gibbs sampler is a Markov chain Monte Carlo (MCMC) algorithm for obtaining a sequence of observations which are approximated from a specified multivariate probability distribution, when direct sampling is difficult. This sequence can be used to approximate the joint distribution (e.g., to generate a histogram of the distribution); to approximate the marginal distribution of one of the variables, or some subset of the variables (for example, the unknown parameters or latent variables); or to compute an integral (such as the expected value of one of the variables). Typically, some of the variables correspond to observations whose values are known, and hence do not need to be sampled.

　　在统计学上，吉布斯取样或吉布斯取样器是马尔科夫链蒙特卡罗（MCMC）算法，用于获得从特定的多变量 概率分布近似的观测序列，当直接取样困难时。该序列可用于近似联合分布（例如，生成分布的直方图）; 近似一个变量的边际分布，或变量的一些子集（例如，未知参数或潜在变量）; 或以计算积分（如预期值的变量之一的）。通常，一些变量对应于其值是已知的观察值，因此不需要进行采样。Gibbs抽样是常用的一种手段统计推断，特别是贝叶斯推理。它是一种随机算法（即使用随机数的算法），并且是用于统计推断的确定性算法的替代，例如期望最大化算法（EM）。如同其它MCMC算法，Gibbs抽样产生马尔可夫链样品，其中的每一个的相关性与附近的样品。因此，如果需要独立的样品，则必须小心（通常通过仅取每个第n个值，例如每100个值）稀释所得到的样品链。此外，从链的开始（老化期）开始的样本可能不能准确地表示所需的分布。

 

　　MCMC是用于构建 Markov chain随机概率分布的抽样的一类算法。MCMC有很多算法，其中比较流行的是Metropolis-Hastings Algorithm，Gibbs Sampling是Metropolis-Hastings Algorithm的一种特殊情况。
　　Markov chain 是一组事件的集合，在这个集合中，事件是一个接一个发生的，并且下一个事件的发生，只由当前发生的事件决定。用数学符号表示就是：
　　　　A={ a1,a2 … ai, ai+1,… at }
　　　　P(ai+1| a1,a2,…ai) = P(ai+1| ai)
　　这里的ai不一定是一个数字，它有可能是一个向量，或者一个矩阵，例如我们比较感兴趣的问题里ai=（g, u, b）这里g表示基因的效应，u表示环境效应，b表示固定效应，假设我们研究的一个群体，g，u，b的联合分布用π（a）表示。事实上，我们研究QTL，就是要找到π（a），但是有时候π（a）并不是那么好找的，特别是我们要估计的a的参数的个数多于研究的个体数的时候。用一般的least square往往效果不是那么好。
解决方案：
　　用一种叫Markov chain Monte Carlo (MCMC)的方法产生Markov chain，产生的Markov chain{a1,a2 … ai, ai+1,… at }具有如下性质：当t 很大时，比如10000，那么at ~ π（a），这样的话如果我们产生一个markov chain：{a1,a2 … ai, ai+1,… a10000}，那么我们取后面9000个样本的平均 a_hat = (g_hat,u_hat,b_hat) = ∑ai / 9000 (i=1001,1002, … 10000)这里g_hat, u_hat, b_hat 就是基因效应，环境效应，以及固定效应的估计值
　　MCMC算法的关键是两个函数：
　　1） q（ai, ai+1），这个函数决定怎么基于ai得到ai+1
　　2） α（ai, ai+1），这个函数决定得到的ai+1是否保留
　　目的是使得at的分布收敛于π（a）
　　Gibbs Sampling的算法：
　　一般来说我们通常不知道π（a），但我们可以得到p（g | u , b）,p(u | g , b), p ( b | g, u )即三个变量的posterior distribution
　　Step1: 给g, u, b 赋初始值：（g0,u0,b0）
　　Step2: 利用p (g | u0, b0) 产生g1
　　Step3: 利用p (u | g1, b0) 产生u1
　　Step4: 利用p (b | g1, u1) 产生b1
　　Step5: 重复step2~step5 这样我们就可以得到一个markov chain {a1,a2 … ai, ai+1,… at}
　　这里的q（ai, ai+1）= p（g | u , b）* p(u | g , b)* p ( b | g, u )

　　

　　For Example:

　　这里通俗点的解释一下。首先，什么是sampling.sampling就是以一定的概率分布,看发生什么事件.举一个例子.甲只能E：吃饭、学习、打球,时间T：上午、下午、晚上,天气W：晴朗、刮风、下雨。现在要一个sample,这个sample可以是：打球+下午+晴朗..

　　问题是我们不知道p(E,T,W),或者说,不知道三件事的联合分布.当然,如果知道的话,就没有必要用gibbs sampling了.但是,我们知道三件事的conditional distribution.也就是说,p(E|T,W),p(T|E,W),p(W|E,T).现在要做的就是通过这三个已知的条件分布,再用gibbs sampling的方法,得到joint distribution.

　　具体方法.首先随便初始化一个组合,i.e.学习+晚上+刮风,然后依条件概率改变其中的一个变量.具体说,假设我们知道晚上+刮风,我们给E生成一个变量,比如,学习-》吃饭.我们再依条件概率改下一个变量,根据学习+刮风,把晚上变成上午.类似地,把刮风变成刮风（当然可以变成相同的变量）.这样学习+晚上+刮风-》吃饭+上午+刮风.

　　同样的方法,得到一个序列,每个单元包含三个变量,也就是一个马尔可夫链.然后跳过初始的一定数量的单元（比如100个）,然后隔一定的数量取一个单元（比如隔20个取1个）.这样sample到的单元,是逼近联合分布的.