MCMC：Gibbs 采样（matlab 实现）

MCMC: The Gibbs Sampler
多元高斯分布的边缘概率和条件概率
Marginal and conditional distributions of multivariate normal distribution

clear, clc
rng('default')

num_samples = 5000;
num_dims = 2;

mu = [0, 0];
rho(1) = .8; rho(2) = .8;

prop_sigma = 1;
minn = [-3, -3]; maxx = [3, 3];
x = zeros(num_samples, num_dims);

x(1, 1) = unifrnd(minn(1), maxx(1));
x(1, 2) = unifrnd(minn(2), maxx(2));

t = 1;
dims = 1:num_dims;
while t < num_samples
    t = t + 1;
    T = [t-1, t];           % 时刻信息的维护，T(1)：上一时刻，T(2)：下一时刻
    for iD = 1:num_dims
        not_idx = (dims ~= iD);
        mu_cond = mu(iD) + rho(iD)*(x(T(iD), not_idx) - mu(not_idx));
        sigma_cond = sqrt(1-rho(iD)^2);
        x(t, iD) = normrnd(mu_cond, sigma_cond);
    end
end

figure;
h1 = scatter(x(:, 1), x(:, 2), 'r.');
hold on

for t=1:50
    plot([x(t, 1), x(t+1, 1)], [x(t, 2), x(t, 2)], 'k-');           % x 轴方向移动，
    plot([x(t+1, 1), x(t+1, 1)], [x(t, 2), x(t+1, 2)], 'k-');       % y 轴方向移动；
    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