gibbs抽样方法学习笔记

　　Gibbs抽样方法是 Markov Chain Monte Carlo（MCMC）方法的一种，也是应用最为广泛的一种。wikipedia称gibbs抽样为

　　In statistics and in statistical physics, Gibbs sampling or a Gibbs sampler is a Markov chain Monte Carlo (MCMC) algorithm for obtaining a sequence of observations which are approximately from a specified multivariate probability distribution (i.e. from the joint probability distribution of two or more random variables), when direct sampling is difficult.

　　意思是，在统计学和统计物理学中，gibbs抽样是马尔可夫链蒙特卡尔理论（MCMC）中用来获取一系列近似等于指定多维概率分布（比如2个或者多个随即变量的联合概率分布）观察样本的算法。

　　MCMC是用于构建 Markov chain随机概率分布的抽样的一类算法。MCMC有很多算法，其中比较流行的是Metropolis-Hastings Algorithm，Gibbs Sampling是Metropolis-Hastings Algorithm的一种特殊情况。

　　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算法的关键是两个函数：

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 )