当前位置：首页 > news >正文

利用MCMC方法产生平稳的马尔科夫链

news 2025/9/22 9:43:05

马尔科夫链蒙特卡洛(MCMC)方法是一种通过构建平稳马尔科夫链来从复杂概率分布中采样的强大技术。

% 利用MCMC方法产生平稳的马尔科夫链
clear; clc; close all;%% 1. 参数设置
fprintf('设置MCMC参数...\n');% MCMC参数
n_samples = 10000;      % 总样本数
burn_in = 2000;         % 预热期样本数
thin = 5;               % 稀释间隔（每thin个样本取1个）% 目标分布参数（多元正态分布）
mu = [1, 2];            % 均值向量
sigma = [1, 0.8; 0.8, 1]; % 协方差矩阵% 建议分布参数（正态分布）
proposal_std = [0.5, 0.5]; % 建议分布的标准差%% 2. 初始化马尔科夫链
fprintf('初始化马尔科夫链...\n');% 初始状态
current_state = [0, 0]; % 初始点% 存储链
chain = zeros(n_samples, length(mu));
acceptance = zeros(n_samples, 1); % 记录是否接受%% 3. 定义目标分布和建议分布
fprintf('定义概率分布函数...\n');% 目标分布（多元正态分布）
target_pdf = @(x) mvnpdf(x, mu, sigma);% 建议分布（对称的，因此Metropolis-Hastings比率为1）
proposal_pdf = @(x, y) prod(normpdf(x, y, proposal_std));% 建议分布采样函数
proposal_rnd = @(y) y + normrnd(0, proposal_std);%% 4. Metropolis-Hastings算法
fprintf('运行Metropolis-Hastings算法...\n');for i = 1:n_samples% 从建议分布中生成候选点candidate = proposal_rnd(current_state);% 计算接受概率acceptance_ratio = target_pdf(candidate) / target_pdf(current_state);% 由于建议分布是对称的，q(candidate|current)/q(current|candidate)=1% 对于非对称建议分布，需要乘以建议分布比率% 决定是否接受候选点if rand() < min(1, acceptance_ratio)current_state = candidate;acceptance(i) = 1;end% 存储当前状态chain(i, :) = current_state;% 显示进度if mod(i, 1000) == 0fprintf('已完成 %d/%d 次迭代\n', i, n_samples);end
end% 计算接受率
acceptance_rate = mean(acceptance) * 100;
fprintf('接受率: %.2f%%\n', acceptance_rate);%% 5. 后处理：去除预热期和稀释
fprintf('处理后处理...\n');% 去除预热期
chain = chain(burn_in+1:end, :);
acceptance = acceptance(burn_in+1:end);% 稀释
chain = chain(1:thin:end, :);
acceptance = acceptance(1:thin:end);fprintf('处理后样本数: %d\n', size(chain, 1));%% 6. 收敛诊断
fprintf('进行收敛诊断...\n');% 绘制轨迹图
figure;
subplot(2, 2, 1);
plot(chain(:, 1));
xlabel('迭代次数');
ylabel('x_1');
title('马尔科夫链轨迹 (x_1)');
grid on;subplot(2, 2, 2);
plot(chain(:, 2));
xlabel('迭代次数');
ylabel('x_2');
title('马尔科夫链轨迹 (x_2)');
grid on;% 绘制自相关函数
subplot(2, 2, 3);
autocorr(chain(:, 1), 50);
title('x_1的自相关函数');subplot(2, 2, 4);
autocorr(chain(:, 2), 50);
title('x_2的自相关函数');% 计算Gelman-Rubin统计量（需要多条链）
fprintf('计算Gelman-Rubin统计量...\n');
n_chains = 4; % 多条链的数量
gr_stats = gelman_rubin(chain, n_chains, mu, sigma, proposal_std);
fprintf('Gelman-Rubin统计量: x_1=%.4f, x_2=%.4f\n', gr_stats(1), gr_stats(2));%% 7. 评估采样质量
fprintf('评估采样质量...\n');% 绘制样本分布与目标分布的比较
figure;% 创建网格用于绘制目标分布
x1 = linspace(min(chain(:,1))-1, max(chain(:,1))+1, 100);
x2 = linspace(min(chain(:,2))-1, max(chain(:,2))+1, 100);
[X1, X2] = meshgrid(x1, x2);
X = [X1(:), X2(:)];% 计算目标分布
Z = mvnpdf(X, mu, sigma);
Z = reshape(Z, size(X1));% 绘制等高线图
subplot(2, 2, 1);
contour(X1, X2, Z, 20);
hold on;
plot(chain(:,1), chain(:,2), '.', 'MarkerSize', 2, 'Color', [0.7, 0.7, 0.7]);
xlabel('x_1');
ylabel('x_2');
title('样本与目标分布');
colorbar;% 绘制边际分布
subplot(2, 2, 2);
histogram(chain(:,1), 50, 'Normalization', 'pdf');
hold on;
plot(x1, normpdf(x1, mu(1), sqrt(sigma(1,1))), 'r', 'LineWidth', 2);
xlabel('x_1');
ylabel('密度');
title('x_1的边际分布');
legend('样本', '真实分布');subplot(2, 2, 3);
histogram(chain(:,2), 50, 'Normalization', 'pdf');
hold on;
plot(x2, normpdf(x2, mu(2), sqrt(sigma(2,2))), 'r', 'LineWidth', 2);
xlabel('x_2');
ylabel('密度');
title('x_2的边际分布');
legend('样本', '真实分布');% 计算样本均值和协方差
sample_mean = mean(chain);
sample_cov = cov(chain);fprintf('真实均值: [%.4f, %.4f]\n', mu(1), mu(2));
fprintf('样本均值: [%.4f, %.4f]\n', sample_mean(1), sample_mean(2));
fprintf('真实协方差矩阵: [%.4f, %.4f; %.4f, %.4f]\n', ...sigma(1,1), sigma(1,2), sigma(2,1), sigma(2,2));
fprintf('样本协方差矩阵: [%.4f, %.4f; %.4f, %.4f]\n', ...sample_cov(1,1), sample_cov(1,2), sample_cov(2,1), sample_cov(2,2));%% 8. 辅助函数定义function gr_stats = gelman_rubin(chain, n_chains, mu, sigma, proposal_std)% 计算Gelman-Rubin统计量% 将现有链分割成多个部分，模拟多条链n_samples = size(chain, 1);chain_length = floor(n_samples / n_chains);chains = cell(n_chains, 1);% 生成多条链for i = 1:n_chains% 使用不同的初始值initial_state = normrnd(0, 1, 1, 2);% 运行MCMCchain_i = mcmc_chain(chain_length, initial_state, mu, sigma, proposal_std);chains{i} = chain_i;end% 计算Gelman-Rubin统计量% 对于每个参数gr_stats = zeros(1, 2);for param = 1:2% 计算链内方差W = 0;for i = 1:n_chainsW = W + var(chains{i}(:, param));endW = W / n_chains;% 计算链均值chain_means = zeros(n_chains, 1);for i = 1:n_chainschain_means(i) = mean(chains{i}(:, param));end% 计算链间方差B = chain_length * var(chain_means);% 计算估计方差var_estimate = (chain_length - 1)/chain_length * W + 1/chain_length * B;% 计算Gelman-Rubin统计量gr_stats(param) = sqrt(var_estimate / W);end
endfunction chain = mcmc_chain(n_samples, initial_state, mu, sigma, proposal_std)% 运行MCMC生成一条链target_pdf = @(x) mvnpdf(x, mu, sigma);proposal_rnd = @(y) y + normrnd(0, proposal_std);current_state = initial_state;chain = zeros(n_samples, length(mu));for i = 1:n_samples% 从建议分布中生成候选点candidate = proposal_rnd(current_state);% 计算接受概率acceptance_ratio = target_pdf(candidate) / target_pdf(current_state);% 决定是否接受候选点if rand() < min(1, acceptance_ratio)current_state = candidate;end% 存储当前状态chain(i, :) = current_state;end
end%% 9. 自适应MCMC算法（可选）
fprintf('尝试自适应MCMC算法...\n');% 自适应MCMC参数
n_adaptive = 5000;          % 自适应阶段样本数
initial_std = [1.0, 1.0];   % 初始建议分布标准差
target_acceptance = 0.234;  % 目标接受率（最优值）% 初始化
adaptive_chain = zeros(n_adaptive, length(mu));
adaptive_acceptance = zeros(n_adaptive, 1);
current_std = initial_std;
current_state = [0, 0];% 运行自适应MCMC
for i = 1:n_adaptive% 从建议分布中生成候选点candidate = current_state + normrnd(0, current_std);% 计算接受概率acceptance_ratio = target_pdf(candidate) / target_pdf(current_state);% 决定是否接受候选点if rand() < min(1, acceptance_ratio)current_state = candidate;adaptive_acceptance(i) = 1;end% 存储当前状态adaptive_chain(i, :) = current_state;% 自适应调整建议分布标准差if mod(i, 100) == 0current_acceptance = mean(adaptive_acceptance(i-99:i));% 根据接受率调整标准差if current_acceptance > target_acceptancecurrent_std = current_std * 1.1; % 增加探索elsecurrent_std = current_std * 0.9; % 减少探索endend
end% 绘制自适应过程
figure;
subplot(2, 1, 1);
plot(adaptive_chain(:,1));
xlabel('迭代次数');
ylabel('x_1');
title('自适应MCMC轨迹 (x_1)');
grid on;subplot(2, 1, 2);
plot(adaptive_chain(:,2));
xlabel('迭代次数');
ylabel('x_2');
title('自适应MCMC轨迹 (x_2)');
grid on;%% 10. 不同建议分布的比较（可选）
fprintf('比较不同建议分布...\n');% 定义不同的建议分布标准差
proposal_stds = {[0.1, 0.1], [0.5, 0.5], [2.0, 2.0]};
n_comparison = 3000;
colors = ['r', 'g', 'b'];figure;
hold on;for i = 1:length(proposal_stds)% 运行MCMCchain_comp = mcmc_chain(n_comparison, [0, 0], mu, sigma, proposal_stds{i});% 计算接受率acceptance_ratio = sum(diff(chain_comp(:,1)) ~= 0) / (n_comparison - 1);% 绘制轨迹plot(chain_comp(:,1), chain_comp(:,2), '.', 'Color', colors(i), 'MarkerSize', 4);fprintf('建议分布标准差 [%.1f, %.1f] - 接受率: %.2f%%\n', ...proposal_stds{i}(1), proposal_stds{i}(2), acceptance_ratio*100);
end% 绘制目标分布等高线
contour(X1, X2, Z, 10, 'k', 'LineWidth', 1.5);
xlabel('x_1');
ylabel('x_2');
title('不同建议分布的比较');
legend('std=0.1', 'std=0.5', 'std=2.0', '目标分布');
grid on;
hold off;%% 11. 输出总结
fprintf('\n===== MCMC模拟总结 =====\n');
fprintf('目标分布: 二元正态分布, μ=[%.2f, %.2f], Σ=[%.2f, %.2f; %.2f, %.2f]\n', ...mu(1), mu(2), sigma(1,1), sigma(1,2), sigma(2,1), sigma(2,2));
fprintf('建议分布: 正态分布, std=[%.2f, %.2f]\n', proposal_std(1), proposal_std(2));
fprintf('总样本数: %d (预热期: %d, 稀释间隔: %d)\n', n_samples, burn_in, thin);
fprintf('最终样本数: %d\n', size(chain, 1));
fprintf('接受率: %.2f%%\n', acceptance_rate);
fprintf('Gelman-Rubin统计量: x_1=%.4f, x_2=%.4f\n', gr_stats(1), gr_stats(2));
fprintf('样本均值: [%.4f, %.4f]\n', sample_mean(1), sample_mean(2));
fprintf('样本协方差: [%.4f, %.4f; %.4f, %.4f]\n', ...sample_cov(1,1), sample_cov(1,2), sample_cov(2,1), sample_cov(2,2));