MH方法生成MCMC链多各不相关变量MATLAB例子
时间: 2023-07-16 18:12:09 浏览: 41
这里给您提供一个简单的 MATLAB 例子来生成多个不相关的 MCMC 链,使用的是 Metropolis-Hastings (MH) 算法。
我们首先生成一个简单的目标分布,例如二元高斯分布。然后我们使用 MH 算法来采样这个分布并生成多个不相关的 MCMC 链。在这个例子中,我们将生成 3 个链,每个链的长度为 10,000。
以下是代码实现:
```matlab
% 目标分布(二元高斯分布)
mu = [1 2];
sigma = [1 0.5; 0.5 2];
target_dist = @(x) mvnpdf(x, mu, sigma);
% MH 算法参数
num_chains = 3; % 链的数量
chain_length = 10000; % 链的长度
proposal_std = 0.5; % 提议分布的标准差
% 初始化 MCMC 链
chains = zeros(chain_length, 2, num_chains);
accept_counts = zeros(num_chains, 1);
% 对于每个链进行采样
for chain_num = 1:num_chains
% 初始化链的起始点
current_position = [0 0];
current_likelihood = target_dist(current_position);
chains(1,:,chain_num) = current_position;
% 进行 MCMC 链采样
for i = 2:chain_length
% 从提议分布中生成一个新的样本
proposal = current_position + randn(1,2)*proposal_std;
% 计算接受率
proposal_likelihood = target_dist(proposal);
acceptance_prob = min(1, proposal_likelihood/current_likelihood);
% 决定是否接受新样本
if rand < acceptance_prob
current_position = proposal;
current_likelihood = proposal_likelihood;
accept_counts(chain_num) = accept_counts(chain_num) + 1;
end
% 将当前样本添加到链中
chains(i,:,chain_num) = current_position;
end
end
% 输出每个链的接受率
disp(['Acceptance rates:', num2str(accept_counts/chain_length)]);
```
这个代码将生成 3 个不相关的 MCMC 链并输出每个链的接受率。您可以使用这些链来估计目标分布的统计量,例如均值和方差。