MCMC：Gibbs 采样（matlab 实现）

好的，我可以为您提供一些关于 Gibbs 采样的介绍和 Matlab 实现的代码。

Gibbs 采样是一种马尔可夫蒙特卡罗方法，用于从联合分布中采样。它是一种基于条件分布的采样方法，每次采样时只更新一个变量，且更新后的样本仍然满足联合分布。

以下是 Matlab 实现的代码示例：

假设我们要从以下联合分布中采样：

p(x,y) = 1/Z * exp(-x^2/2 - y^2/2 - xy/3)

其中，Z 是归一化常数。

首先，我们可以定义一个初始化的样本点：

x_init = 0;
y_init = 0;

然后，我们可以设置采样次数和 Burn-in 阶段的样本数：

n_samples = 10000; % 采样次数
n_burnin = 1000; % Burn-in 阶段的样本数

接下来，我们可以编写 Gibbs 采样的代码：

samples = zeros(n_samples, 2);
x = x_init;
y = y_init;

for i = 1:n_samples+n_burnin
% 从条件分布 p(x|y) 中采样
x = normrnd(y/3, sqrt(1/2));

% 从条件分布 p(y|x) 中采样
y = normrnd(x/3, sqrt(1/2));

if i > n_burnin
samples(i-n_burnin, :) = [x, y];
end
end

最后，我们可以绘制样本点的散点图和联合分布的等高线图：

figure;
scatter(samples(:, 1), samples(:, 2), 10, '.');
xlabel('x');
ylabel('y');
title('Scatter Plot');

figure;
x_min = min(samples(:, 1));
x_max = max(samples(:, 1));
y_min = min(samples(:, 2));
y_max = max(samples(:, 2));
x_range = linspace(x_min, x_max, 100);
y_range = linspace(y_min, y_max, 100);
[X,Y] = meshgrid(x_range, y_range);
Z = exp(-X.^2/2 - Y.^2/2 - X.*Y/3);
contour(X, Y, Z);
xlabel('x');
ylabel('y');
title('Contour Plot');

希望这个示例代码能够帮助您理解 Gibbs 采样的实现过程。