吉布斯采样——原理及matlab实现

Posted 2022-08-13 风翼冰舟

原文来自：https://victorfang.wordpress.com/2014/04/29/mcmc-the-gibbs-sampler-simple-example-w-matlab-code/

【注】评论区有同学指出译文理论编码有误，请参考更官方的文献，个人当时仅验证过红色字体部分理论与维基百科中二位随机变量吉布斯采样的结果是否对应，其余部分有意见希望可以详细指出，大家互相交流。

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个符合分布的点

 

%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