中国站长工具,展览展示设计必看网站,有网站是做水果原产地代发的吗,建设银行网站怎么设置转账额度原文来自#xff1a;https://victorfang.wordpress.com/2014/04/29/mcmc-the-gibbs-sampler-simple-example-w-matlab-code/
【注】评论区有同学指出译文理论编码有误#xff0c;请参考更官方的文献#xff0c;个人当时仅验证过红色字体部分理论与维基百科中二位随机变量吉…原文来自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.设置t0
2.生成初始状态
3.重复直到tM { 对于每一个维度i1...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 SAMPLESnIx dims~iD; % *NOT* THE CURRENT DIMENSION找到另外一维nIx[0 1]logical类型% CONDITIONAL MEANmuCond mu(iD) rho(iD)*(x(T(iD),nIx)-mu(nIx));%计算均值表达式μ(1)ρ(1)*(x(n,2)-μ(2))其中x(n,2)代表样本第n个数据的第二维% CONDITIONAL VARIANCEvarCond sqrt(1-rho(iD)^2);%计算方差% DRAW FROM CONDITIONALx(t,iD) normrnd(muCond,varCond);%正态分布随机函数计算得到当前第t个数据的第1维数据valueend
end% DISPLAY SAMPLING DYNAMICS
figure;
h1 scatter(x(:,1),x(:,2),r.);%scatter描绘散点图x为横坐标y为纵坐标% CONDITIONAL STEPS/SAMPLES
hold on;%画出前五十个采样点
for t 1:50plot([x(t,1),x(t1,1)],[x(t,2),x(t,2)],k-);plot([x(t1,1),x(t1,1)],[x(t,2),x(t1,2)],k-);h2 plot(x(t1,1),x(t1,2),ko);
endh3 scatter(x(1,1),x(1,2),go,Linewidth,3);
legend([h1,h2,h3],{Samples,1st 50 Samples,x(t0)},Location,Northwest)
hold off;
xlabel(x_1);
ylabel(x_2);
axis square
