MH方法生成MCMC链多变量MATLAB例子

MH方法是基于Metropolis-Hastings算法的一种MCMC采样方法，用于生成多变量随机分布的样本。下面是一个MATLAB例子，用于演示如何使用MH方法生成多变量样本。

假设我们要生成一个二元正态分布的样本，可以按照以下步骤进行：

1. 定义目标分布函数

% 目标分布函数
mu = [0, 0];
Sigma = [1, 0.5; 0.5, 1];
mvnpdf = @(x) mvnpdf(x, mu, Sigma);

其中，mu是分布的均值，Sigma是协方差矩阵，mvnpdf是MATLAB自带的多元正态分布概率密度函数。

2. 初始化样本

% 初始化样本
x = randn(1, 2);

我们随机生成一个初始样本，其服从标准正态分布。

3. 定义提议分布函数

% 提议分布函数
proposal = @(x) mvnrnd(x, 0.5*eye(2));

我们选择正态分布作为提议分布函数，其均值为当前样本，方差为0.5。

4. 生成样本

% 生成样本
N = 10000;
samples = zeros(N, 2);

for i = 2:N
    % 生成提议样本
    x_star = proposal(x(i-1, :));
    
    % 计算接受概率
    alpha = min(1, mvnpdf(x_star)/mvnpdf(x(i-1, :))*...
        normpdf(x(i-1, :), x_star, 0.5*eye(2))/...
        normpdf(x_star, x(i-1, :), 0.5*eye(2)));
    
    % 接受或拒绝提议样本
    if rand < alpha
        x(i, :) = x_star;
    else
        x(i, :) = x(i-1, :);
    end
end

samples = x;

我们根据MH算法，迭代生成10000个样本，其中每个样本都是从当前样本中生成的提议样本中接受或拒绝的。最终得到的样本即为目标分布的样本。

5. 绘制结果

% 绘制结果
figure;
scatter(samples(:, 1), samples(:, 2), '.');
title('MH Sampled Distribution');
xlabel('x1');
ylabel('x2');

最终得到的结果如下图所示：