20171208更新:
之前理解有错,没有体现出Gibbs抽样的“抽样”操作。重贴主程序代码如下:
close all; clear all; rng(0);
%% Gibbs sampling method for Ising model
% chapter 24.2.2 in MLaPP
% by forseerwang, http://blog.csdn.net/foreseerwang
% QQ: 50834
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Image generation
a=double(imread('LetterB.bmp'))*2-1;
figure(1); subplot(2,2,1);
imagesc(a);
colorbar;
title('original image', 'fontsize', 20);
axis('square');
colormap gray;
axis off;
% add noise
sigma = 2; % noise level
y = a + sigma*randn(size(a)); %y = noisy signal
figure(1); subplot(2,2,2);
imagesc(y);
colorbar;
title('noisy image', 'fontsize', 20);
axis('square');
colormap gray;
axis off;
%%
% 迭代参数
maxIter = 30;
% 最初需要舍弃的抽样数
Nburnin = 10;
psi_p = log(normpdf(y(:), 1, sigma));
psi_m = log(normpdf(y(:), -1, sigma));
logodds = psi_p - psi_m;
J = 1;
xpred = y(:);
xrec = zeros(length(xpred), maxIter - Nburnin);
[M,N] = size(y);
for iter = 1:maxIter
for ii = 1:N,
for jj = 1:M,
pos = jj + M*(ii-1);
nbr = pos + [-1,1,-M,M];
nbr(([jj==1,jj==M,ii==1,ii==N])) = [];
eta = sum(xpred(nbr));
% MLaPP中的公式(24.6),书中有typo
prob = fun_sigm(2*J*eta + logodds(pos));
% 抽样
if rand>prob,
xpred(pos) = -1;
else
xpred(pos) = +1;
end;
end;
end;
if iter>Nburnin,
xrec(:,iter-Nburnin) = xpred;
end;
end;
xmean = mean(xrec,2);
figure(1); subplot(2,2,3);
imagesc(reshape(xmean,M,N));
colorbar;
title(sprintf('Mean image over last %d iter', maxIter-Nburnin), 'fontsize', 20);
axis('square');
colormap gray;
axis off;
xpred=reshape(xpred,M,N);
figure(1); subplot(2,2,4);
imagesc(xpred);
colorbar;
title(sprintf('image after %d iter', maxIter), 'fontsize', 20);
axis('square');
colormap gray;
axis off;
这个代码的运行结果如下:
%% 旧的版本,理解和代码有误。
Gibbs SamplingGibbs Sampling是MCMC(蒙特卡洛马尔科夫链:Markov chain Monte Carlo)中最流行的一个算法。抛开数学推导不说,这个算法的形式非常简单。详见MLaPP第24.2节。
作为一个入门示例,此处给出了Gibbs sampling用于Ising模型的实现,即MLaPP第24.2.2节的实现。
算法很简单,核心公式就是(24.1)和(24.6),只是需要注意的是公式(24.6)里有typo,括号中的减号应该是加号,这个自己稍微推导一下就出来了,很简单。
还是我之前画的那个“B”,降噪结果如下(可以和之前Mean Field算法对比一下,http://blog.csdn.net/foreseerwang/article/details/78343238):
%% 20171208更新:以下主程序代码有误。
代码如下:
close all; clear all; rng(0);
%% Gibbs sampling method for Ising model
% chapter 24.2.2 in MLaPP
% by forseerwang, http://blog.csdn.net/foreseerwang
% QQ: 50834
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Image generation
a=double(imread('LetterB.bmp'))*2-1;
figure(1); subplot(2,2,1);
imagesc(a);
colorbar;
title('original image', 'fontsize', 20);
axis('square');
colormap gray;
axis off;
% add noise
sigma = 2; % noise level
y = a + sigma*randn(size(a)); %y = noisy signal
figure(1); subplot(2,2,2);
imagesc(y);
colorbar;
title('noisy image', 'fontsize', 20);
axis('square');
colormap gray;
axis off;
%%
% 迭代参数
maxIter = 10;
psi_p = log(normpdf(y(:), 1, sigma));
psi_m = log(normpdf(y(:), -1, sigma));
logodds = psi_p - psi_m;
J = 1;
xpred = y(:);
[M,N] = size(y);
for iter = 1:maxIter
for ii = 1:N,
for jj = 1:M,
pos = jj + M*(ii-1);
nbr = pos + [-1,1,-M,M];
nbr(([jj==1,jj==M,ii==1,ii==N])) = [];
eta = sum(xpred(nbr));
% MLaPP中的公式(24.6),书中有typo
xpred(pos) = fun_sigm(2*J*eta + logodds(pos))*2-1;
end;
end;
if iter == round(maxIter/2),
xpred=reshape(xpred,M,N);
figure(1); subplot(2,2,3);
imagesc(xpred);
colorbar;
title(sprintf('image after %d iter', iter), 'fontsize', 20);
axis('square');
colormap gray;
axis off;
end;
end;
xpred=reshape(xpred,M,N);
figure(1); subplot(2,2,4);
imagesc(xpred);
colorbar;
title(sprintf('image after %d iter', maxIter), 'fontsize', 20);
axis('square');
colormap gray;
axis off;