Gibbs sampling

关于Gibbs sampling，首先看一下Wiki上的解释：Gibbs sampling or Gibbs sampler is an algorithm to generate a sequence of samples from the joint probability distribution of two or more random variables. The purpose of such a sequence is to approximate the joint distribution, or to compute an integral (such as an expected value).

说到Gibbs Sampling 就不得不说markov chain了。

Markov chain 是一组事件的集合，在这个集合中，事件是一个接一个发生的，并且下一个事件的发生，只由当前发生的事件决定。用数学符号表示就是：

A={ a₁,a₂ … a_i, a_i+1,… a_t }

P(a_i+1| a₁,a₂,…a_i) = P(a_i+1| a_i)

这里的a_i不一定是一个数字，它有可能是一个向量，或者一个矩阵，例如我们比较感兴趣的问题里a_i=（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{a₁,a₂ … a_i, a_i+1,… a_t }具有如下性质：当t 很大时，比如10000，那么a_t ~ π（a），这样的话如果我们产生一个markov chain：{a₁,a₂ … a_i, a_i+1,… a₁₀₀₀₀}，那么我们取后面9000个样本的平均

a_hat = (g_hat,u_hat,b_hat) = ∑a_i/ 9000 (i=1001,1002, … 10000)

这里g_hat, u_hat, b_hat 就是基因效应，环境效应，以及固定效应的估计值

MCMC有很多算法，其中比较流行的是Metropolis-Hastings Algorithm，Gibbs Sampling是Metropolis-Hastings Algorithm的一种特殊情况。MCMC算法的关键是两个函数：

1） q（a_i, a_i+1），这个函数决定怎么基于a_i得到a_i+1

2） α（a_i, a_i+1），这个函数决定得到的a_i+1是否保留

目的是使得a_t的分布收敛于π（a）

Gibbs Sampling的算法：

一般来说我们通常不知道π（a），但我们可以得到p（g | u , b）,p(u | g , b), p ( b | g, u )即三个变量的posterior distribution

Step1: 给g, u, b 赋初始值：（g₀,u₀,b₀）

Step2: 利用p (g | u₀, b₀) 产生g1

Step3: 利用p (u | g₁, b₀) 产生u1

Step4: 利用p (b | g₁, u₁) 产生b1

Step5: 重复step2~step5 这样我们就可以得到一个markov chain {a₁,a₂ … a_i, a_i+1,… a_t}

这里的q（a_i, a_i+1）= p（g | u , b）* p(u | g , b)* p ( b | g, u )

这里通俗点的解释一下。首先，什么是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到的单元，是逼近联合分布的。

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

互动百科这么说的：

概率推理的通用方法，是Metropolis-Hastings算法的一个特例，因此也是Markov chain Monte Carlo算法的一种。

虽然它的通用性比较好，但导致了计算代价较高，所以在许多应用里，包括具有不完备信息的应用，都采用其它更为高效的方法。然而，理解这一方法有助于增进对统计推理问题的理解。

中心思想

由一个具有2个或更多变量的联合概率分布P(x1,x2,…,xn)，生成一个样本序列{y1,y2,…,ym}，用于逼近这一个联合分布，或计算一个积分（例如期望）。

适用于处理不完备信息，当联合分布不明确，而各个变量的条件分布已知的情况。

根据其他变量的当前值，依次对分布的每个变量生成一个实例。

随机过程

对一个随机过程，例如马尔可夫链过程，一般包括一个有限的状态集合和一个概率转移矩阵。假设这个过程各个各个状态都是可遍历的(ergodic)，即转移矩阵中的元素值都大于0。为此，我们可以选择任意状态为初始态 Q0，计算转化N次后可能到达的状态 Qn 的概率。当N取值足够大时，可以计算得到这一过程最有可能的终态。

假设有一个变量集合X={X1，X2，……，Xn}，P(X)为集合X的联合分布，0<p(x)<1。< p="">

我们将这些变量看做一个马尔科夫过程中的状态集，这一过程定义为：S=∏i=1~n