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

马尔科夫链蒙特卡洛(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) == 0
        fprintf('已完成 %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);
        
        % 运行MCMC
        chain_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_chains
            W = W + var(chains{i}(:, param));
        end
        W = W / n_chains;
        
        % 计算链均值
        chain_means = zeros(n_chains, 1);
        for i = 1:n_chains
            chain_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
end

function 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) == 0
        current_acceptance = mean(adaptive_acceptance(i-99:i));
        
        % 根据接受率调整标准差
        if current_acceptance > target_acceptance
            current_std = current_std * 1.1; % 增加探索
        else
            current_std = current_std * 0.9; % 减少探索
        end
    end
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)
    % 运行MCMC
    chain_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));

程序说明

这个MATLAB程序展示了如何使用MCMC方法（特别是Metropolis-Hastings算法）产生平稳的马尔科夫链，包括以下主要内容：

1. MCMC参数设置

• 定义样本数量、预热期和稀释间隔
• 设置目标分布（多元正态分布）的参数
• 配置建议分布（正态分布）的参数

2. Metropolis-Hastings算法实现

• 初始化马尔科夫链
• 定义目标分布和建议分布函数
• 实现Metropolis-Hastings算法的核心步骤：
  1. 从建议分布生成候选点
  2. 计算接受概率
  3. 根据接受概率决定是否接受候选点
  4. 更新链状态

3. 后处理与收敛诊断

• 去除预热期样本
• 应用稀释以减少自相关性
• 绘制链轨迹和自相关函数
• 计算Gelman-Rubin统计量评估收敛性

4. 采样质量评估

• 比较样本分布与目标分布
• 计算样本均值和协方差矩阵
• 评估与真实参数的接近程度

5. 高级功能

• 实现自适应MCMC算法，动态调整建议分布参数
• 比较不同建议分布对采样效率的影响
• 提供完整的性能评估和总结报告

参考代码 利用MCMC（马尔科夫蒙特卡洛）方法，产生平稳的马尔科夫链 www.youwenfan.com/contentcnh/64018.html

关键概念

Metropolis-Hastings算法

Metropolis-Hastings算法是MCMC方法中最常用的算法之一，其核心步骤如下：

1. 从建议分布\(q(x'|x_t)\)中抽取候选样本\(x'\)
2. 计算接受概率：\(\alpha = \min\left(1, \frac{p(x')q(x_t|x')}{p(x_t)q(x'|x_t)}\right)\)
3. 以概率\(\alpha\)接受候选样本，设置\(x_{t+1} = x'\)；否则\(x_{t+1} = x_t\)

平稳分布

马尔科夫链的平稳分布是指当链达到稳定状态后，状态的概率分布不再随时间变化。MCMC方法的核心思想是构建一个马尔科夫链，使其平稳分布等于目标分布。

收敛诊断

• 轨迹图：直观检查链是否达到稳定状态
• 自相关函数：评估样本之间的相关性
• Gelman-Rubin统计量：通过比较多条链的方差评估收敛性

使用说明

1. 修改目标分布参数以适应您的具体问题
2. 调整建议分布的标准差以优化接受率（理想接受率约为23.4%）
3. 根据链的混合情况调整预热期和稀释间隔
4. 使用收敛诊断工具确保链已达到平稳状态

扩展应用

这个MCMC框架可以扩展到各种复杂分布和问题：

1. 高维分布：处理更高维度的目标分布
2. 非正态分布：适应各种复杂的概率分布形式
3. 分层贝叶斯模型：应用于复杂的统计模型
4. 模型选择：结合贝叶因子进行模型比较
5. 实时MCMC：实现在线学习和自适应采样

注意事项

1. 建议分布的选择对MCMC效率至关重要
2. 预热期长度应足够长，确保链达到平稳状态
3. 稀释可以减少样本间的自相关性，但会增加计算成本
4. 对于复杂分布，可能需要使用更先进的MCMC算法（如HMC、NUTS等）

这个实现提供了一个完整的MCMC框架，您可以根据实际需求进行修改和扩展，以应用于各种复杂的概率分布采样问题。