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

巴特西

MCMC：Gibbs 采样（matlab 实现）

最新文章

热门文章