吉布斯采样是生成马尔科夫链的一种方法，生成的马尔科夫链可以用来做蒙特卡洛仿真，从而求得一个较复杂的多元分布。
吉布斯采样的具体做法：假设有一个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 是一组事件的集合，在这个集合中，事件是一个接一个发生的，并且下一个事件的发生，只由当前发生的事件决定。用数学符号表示就是：
不一定是一个数字，它有可能是一个向量，或者一个矩阵，例如我们比较感兴趣的问题里
解决方案：
具有如下性质：当t 很大时，比如10000，那么
~ π（a），这样的话如果我们产生一个markov chain：
这里g_hat, u_hat, b_hat 就是基因效应，环境效应，以及固定效应的估计值；
MCMC算法的关键是两个函数：
,
），这个函数决定怎么基于
得到
,
），这个函数决定得到的
的分布收敛于π（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;
注意：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个独立的表达式：
吉布斯采样的工作方法和分成分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实现吉布斯采样，输出的采样如下：
总结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代码
[html] view plain copy
print ?
1. %https://victorfang.wordpress.com/2014/04/29/mcmc-the-gibbs-sampler-simple-example-w-matlab-code/
2. %seed 用来控制 rand 和 randn
3. % 如果没有设置seed，每次运行rand或randn产生的随机数都是不一样的
4. % 用了seed，比如设置rand('seed',0);，那么每次运行rand产生的随机数是一样的，这样对调试程序很有帮助
5. rand('seed' ,12345);
6.
7. nSamples = 5000;
8.
9. mu = [0 0]; % TARGET MEAN目标均值
10. rho(1) = 0.8; % rho_21目标方差
11. rho(2) = 0.8; % rho_12目标方差
12.
13. % INITIALIZE THE GIBBS SAMPLER
14. propSigma = 1; % PROPOSAL VARIANCE
15. minn = [-3 -3];
16. maxx = [3 3];
17.
18. % INITIALIZE SAMPLES
19. x = zeros(nSamples,2);
20. x(1,1) = unifrnd(minn(1), maxx(1));%unifrnd生成连续均匀分布的随机数
21. x(1,2) = unifrnd(minn(2), maxx(2));
22.
23. dims = 1:2; % INDEX INTO EACH DIMENSION
24.
25. % RUN GIBBS SAMPLER
26. t = 1;
27. while t < nSamples%总共采样出5000个采样点
28. t = t + 1;
29. T = [t-1,t];
30. iD = 1:2 % LOOP OVER DIMENSIONS总共两维，注释先讨论第一维
31. % UPDATE SAMPLES
32. nIx = dims~=iD; % *NOT* THE CURRENT DIMENSION找到另外一维nIx=[0 1]logical类型
33. % CONDITIONAL MEAN
34. muCond = mu(iD) + rho(iD)*(x(T(iD),nIx)-mu(nIx));%计算均值=表达式μ(1)+ρ(1)*(x(n,2)-μ(2))其中x(n,2)代表样本第n个数据的第二维
35. % CONDITIONAL VARIANCE
36. varCond = sqrt(1-rho(iD)^2);%计算方差
37. % DRAW FROM CONDITIONAL
38. x(t,iD) = normrnd(muCond,varCond);%正态分布随机函数，计算得到当前第t个数据的第1维数据value
39. end
40. end
41.
42. % DISPLAY SAMPLING DYNAMICS
43. figure;
44. h1 = scatter(x(:,1),x(:,2),'r.');%scatter描绘散点图，x为横坐标，y为纵坐标
45.
46. % CONDITIONAL STEPS/SAMPLES
47. hold on;%画出前五十个采样点
48. for t = 1:50
49. plot([x(t,1),x(t+1,1)],[x(t,2),x(t,2)],'k-');
50. plot([x(t+1,1),x(t+1,1)],[x(t,2),x(t+1,2)],'k-');
51. h2 = plot(x(t+1,1),x(t+1,2),'ko');
52. end
53.
54. h3 = scatter(x(1,1),x(1,2),'go','Linewidth',3);
55. legend([h1,h2,h3],{'Samples','1st 50 Samples','x(t=0)'},'Location','Northwest')
56. hold off;
57. xlabel('x_1');
58. ylabel('x_2');
59. axis square
