吉布斯采样是生成马尔科夫链的一种方法，生成的马尔科夫链可以用来做蒙特卡洛仿真，从而求得一个较复杂的多元分布。

吉布斯采样的具体做法：假设有一个k维的随机向量，现想要构造一条有n个样本的k维向量（n样本马尔科夫序列），那么（随机）初始化一个k维向量，然后固定这个向量其中的k-1个元素，抽取剩下的那个元素（生成给定后验的随机数），这样循环k次，就把整个向量更新了一遍，也就是生成了一个新的样本，把这个整体重复n次就得到了一条马尔科夫链。

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

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

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

这里的

不一定是一个数字，它有可能是一个向量，或者一个矩阵，例如我们比较感兴趣的问题里

=（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

具有如下性质：当t 很大时，比如10000，那么

~ π（a），这样的话如果我们产生一个markov chain：

，那么我们取后面9000个样本的平均：

a_hat = (g_hat,u_hat,b_hat) = ∑

/ 9000 (i=1001,1002, … 10000)；

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

MCMC算法的关键是两个函数：

1）q（

），这个函数决定怎么基于

得到

；

2）α（

），这个函数决定得到的

是否保留；

目的是使得

的分布收敛于π（a） [1] 。

过程

编辑

一般来说我们通常不知道π（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

。

这里的q（

）= p（g | u , b）* p(u | g , b)* p ( b | g, u )

注意：Gibbs采样的目的是获得一个样本，不是计算概率，但可以通过其他方法来统计概率

MCMC(马尔可夫链蒙特卡洛方法)：the Gibbs Sampler(吉布斯采样)

在之前的博客中，我们对比了在一个多元概率分布p(x)中，分别采用分组（block-wise）和分成分（component-wise）实现梅特罗波利斯哈斯廷斯算法。对于多元变量问题中的MCMC算法，分成分更新权重比分组更新更有效，因为通过使每一个成分/维度独立于其他成分/维度，我们将更可能去接受一个提议采样【注，这个proposed sample应该就是前面博客里面提到的转移提议分布】。然而，提议采样仍然可能被拒绝，导致有些多余的计算，因为他们被拒绝了，计算了但是一直未使用。吉布斯采样是另外一种比较受欢迎的MCMC采样技术，提供了避免这种多余计算的方法。就想分成分实现Metropolis Hastings算法，吉布斯仍然使用分成分更新。然而，不像Metropolis Hastings采样，所有的提议采样将被接受，因此不会有多余的计算。

基于两个标准，吉布斯采样使用某些类别的问题。给定一个目标分布p(x)，其中，第一个标准是以其他所有变量联合起来的联合分布为条件的每一个变量的条件分布有解析(数学)表达式。在形式上，如果目标分布p(x)是D维的，我们必须有D个独立的表达式：

这些表达式的每一个都定义了在知道其他j(j≠i)维的数值的情况下第i维的概率。具有每一个变量的条件分布代表我们不需要像Metropolis Hastings算法需要提议分布或者接受/拒绝标准。因此，当其他变量固定的时候，我们可以简单的从每一个条件中去采样。第二个标准就是我们必须能够从每一个条件分布中去采样。如果我们想要去实现一个算法，这个附加条件是非常明显的。

吉布斯采样的工作方法和分成分Metropolis Hastings算法很像，除了取缔借鉴每一个维度的提议分布，然后对于接受或者拒绝提议采样，我们采用简单地依据变量对应的条件分布去选取此维度的值。我们会接受所有选取的值。类似分成分Metropolis Hastings算法，我们依次通过每一个变量，在其它变量固定的时候对它采样。吉布斯采样的步骤大致如下：

1.设置t=0

2.生成初始状态

3.重复直到t=M
{
对于每一个维度i=1...D
从中得到

}
为了对吉布斯采样有更好的理解，我们下面来实现一下吉布斯采样，去解决与前面提到过的同样的多元变量采样问题。

例子：从二元正态分布中采样Example: Sampling from a bivariate a Normal distribution

这个例子与前面一样，从2维的正态分布使用分组和分成分的Metropolis-Hastings算法采样。这里我们展示使用同样的目标分布如何实现吉布斯采样。重复提示一下，目标函数p(x)是一种规范化形式，表示如下：

①均值是

②方差是

为了使用吉布斯采样从这个分布中采样，我们需要有变量/维度x1和x2的条件分布：

是第二个维度的前一个状态，是从中得到的第一个维度的状态。有差异的原因是更新x1和x2用的是(t-1)和t时刻的状态，在上一节中的算法大纲第三步可以看出来。第t次迭代，我们首先以变量x2的最近状态即第(t-1)次迭代结果为条件，为x1采样一种新状态。然后再以第t次迭代得到的x1的最新状态为条件采样得到变量x2。
经过一些数学推导(这里先跳过，下面会有详细的过程)，我们发现目标正态分布的两个条件分布是：

每一个都是单变量的正态分布，其中均值依赖条件变量的最近状态的值，方差依赖两个变量之间的目标方差。
使用上述描述的变量x1和x2的条件概率，我们下面采用matlab实现吉布斯采样，输出的采样如下：

观察上表，发现每一次迭代中吉布斯采样中的马尔可夫链是如何第一次沿着x1方向前进一步，然后沿着x2的方向前进。这展示了吉布斯采样以分成分方法一次单独为每一个变量采样。

总结Wrapping Up

吉布斯采样是为复杂多元概率分布采样的一个受欢迎的MCMC方法。然而，吉布斯采样不能用于一般的抽样问题。对于许多目标分布，很难或者不可能去获取到所有需要的条件分布的近似表达。在其它情况下，对于所有条件的解析表达式或许存在，但是它或许很难从任意的或者全部的条件分布去采样(在这种情况下，使用单变量( univariate sampling methods)采样比如拒绝抽样(rejection sampling)和Metropolis类型的MCMC技术去逼近每一个条件的样本是比较普遍的)。吉布斯采样是非常受欢迎的贝叶斯方法，模型经常以这种方式设计：所有模型变量的条件表达式非常容易获取，并且采用一种能够被高效采样的众所周知的形式。
吉布斯采样，就想很多MCMC方法，有“慢混合(slow mixing)”的问题。慢混合的发生是在潜在的马尔可夫链需要很长时间去充分探索出x的值，从而给出一个更好的p(x)的表征(characterization)。慢混合是因为一些因素包括马尔可夫链的“随机走动(random walk)”特性，并且马尔可夫链有“卡住”的趋势，因为仅仅采样了x的一个单独区域，这个区域在p(x)下有很高的概率。这种反应对于多模式(multiple modes)或者重尾(heavy tails)中的分布进行采样效果不好，比如混合蒙特卡洛已被发展成一个包含附加动力学(incorporate additional dynamics)的能提高马尔可夫链路径效率的方法。将来会讨论混合蒙特卡洛方法

matlab代码

实现的应该是给定了一个均值和方差，以及初始的一个样本点，然后采样出5000个符合分布的点

[html]  view plain copy  
 print ? 
 %https://victorfang.wordpress.com/2014/04/29/mcmc-the-gibbs-sampler-simple-example-w-matlab-code/  
 %seed 用来控制 rand 和 randn   
 % 如果没有设置seed，每次运行rand或randn产生的随机数都是不一样的  
 % 用了seed，比如设置rand('seed',0);，那么每次运行rand产生的随机数是一样的，这样对调试程序很有帮助  
 rand('seed' ,12345);  
   
 nSamples = 5000;  
    
 mu = [0 0]; % TARGET MEAN目标均值  
 rho(1) = 0.8; % rho_21目标方差  
 rho(2) = 0.8; % rho_12目标方差  
    
 % INITIALIZE THE GIBBS SAMPLER  
 propSigma = 1; % PROPOSAL VARIANCE  
 minn = [-3 -3];  
 maxx = [3 3];  
    
 % INITIALIZE SAMPLES  
 x = zeros(nSamples,2);  
 x(1,1) = unifrnd(minn(1), maxx(1));%unifrnd生成连续均匀分布的随机数  
 x(1,2) = unifrnd(minn(2), maxx(2));  
    
 dims = 1:2; % INDEX INTO EACH DIMENSION  
    
 % RUN GIBBS SAMPLER  
 t = 1;  
 while t < nSamples%总共采样出5000个采样点  
     t = t + 1;  
     T = [t-1,t];  
     for iD = 1:2 % LOOP OVER DIMENSIONS总共两维，注释先讨论第一维  
         % UPDATE SAMPLES  
         nIx = dims~=iD; % *NOT* THE CURRENT DIMENSION找到另外一维nIx=[0 1]logical类型  
         % CONDITIONAL MEAN  
         muCond = mu(iD) + rho(iD)*(x(T(iD),nIx)-mu(nIx));%计算均值=表达式μ(1)+ρ(1)*(x(n,2)-μ(2))其中x(n,2)代表样本第n个数据的第二维  
         % CONDITIONAL VARIANCE  
         varCond = sqrt(1-rho(iD)^2);%计算方差  
         % DRAW FROM CONDITIONAL  
         x(t,iD) = normrnd(muCond,varCond);%正态分布随机函数，计算得到当前第t个数据的第1维数据value  
     end  
 end  
    
 % DISPLAY SAMPLING DYNAMICS  
 figure;  
 h1 = scatter(x(:,1),x(:,2),'r.');%scatter描绘散点图，x为横坐标，y为纵坐标  
    
 % CONDITIONAL STEPS/SAMPLES  
 hold on;%画出前五十个采样点  
 for t = 1:50  
     plot([x(t,1),x(t+1,1)],[x(t,2),x(t,2)],'k-');  
     plot([x(t+1,1),x(t+1,1)],[x(t,2),x(t+1,2)],'k-');  
     h2 = plot(x(t+1,1),x(t+1,2),'ko');  
 end  
    
 h3 = scatter(x(1,1),x(1,2),'go','Linewidth',3);  
 legend([h1,h2,h3],{'Samples','1st 50 Samples','x(t=0)'},'Location','Northwest')  
 hold off;  
 xlabel('x_1');  
 ylabel('x_2');  
 axis square  