一、MCMC核心算法实现
1. Metropolis-Hastings算法
MATLAB代码示例（目标分布为双峰混合高斯分布）：
% 参数设置
N = 1e5; % 迭代次数
x = zeros(N,1); % 初始化链
x(1) = 0; % 初始值
mu1 = 1; sigma1 = 1; % 第一个高斯分布参数
mu2 = 5; sigma2 = 1; % 第二个高斯分布参数% 目标分布定义
target = @(x) 0.3*normpdf(x,mu1,sigma1) + 0.7*normpdf(x,mu2,sigma2);% 提议分布（高斯随机游走）
step = 1; % 步长% MCMC迭代
for i = 2:Nx_star = x(i-1) + step*randn; % 生成候选点% 计算接受概率alpha = min(1, target(x_star)/target(x(i-1)));if rand < alphax(i) = x_star;elsex(i) = x(i-1);end
end% 可视化结果
figure;
histogram(x,50,'Normalization','pdf');
hold on;
x_plot = linspace(-2,8,1000);
plot(x_plot,target(x_plot),'r','LineWidth',2);
title('MCMC采样结果与目标分布对比');
xlabel('x'); ylabel('密度');
2. 吉布斯采样（Gibbs Sampling）
应用场景：多元高斯分布参数估计
% 生成模拟数据
mu_true = [2; -1]; Sigma = [1 0.8; 0.8 1];
data = mvnrnd(mu_true,Sigma,1000);% 初始化参数
mu = [0; 0]; sigma = 1;% 吉布斯采样
N = 5000; burn_in = 1000;
for i = 1:N% 更新mu | sigmamu = mvnrnd(mu, inv(sigma)*data'/(inv(sigma)*data'*inv(sigma)+eye(2)/10));% 更新sigma | musigma = wishrnd(inv(data'*data/10 + eye(2)), 3);
end% 去除燃烧期
mu = mu(burn_in:end,:);
二、收敛性诊断方法
1. 图形化诊断工具
% 绘制轨迹图
figure;
plot(x(1:1000:end));
title('马尔可夫链轨迹图');% 自相关分析
lags = 1:50;
ac = autocorr(x,lags);
figure;
stem(lags,ac,'filled');
title('自相关函数');
2. 统计检验指标
% Geweke诊断
[ge,gep](@ref)= geweke(x);
disp(['Geweke检验p值: ', num2str(gep)]);% 有效样本量计算
ess = effective_sample_size(x);
disp(['有效样本量: ', num2str(ess)]);
三、优化
1. 自适应步长调整
% 根据接受率动态调整步长
if acceptance_rate > 0.4step = step * 1.1;
elseif acceptance_rate < 0.2step = step * 0.9;
end
2. 并行计算加速
% 使用parfor实现并行采样
parfor i = 1:num_chainschain(:,:,i) = run_chain(theta_init);
end
3. Hamiltonian Monte Carlo
% 引入动量变量
p = randn(size(theta));
for i = 1:N[p, theta](@ref)= leapfrog(p, theta, epsilon, L);if rand < exp(-H(p,theta))accept(i) = 1;end
end
四、完整代码模板
function mcmc_demo()% 参数初始化theta_init = [0; 0];num_iter = 1e5;burn_in = 1e4;% 构建目标分布target = @(x) log_target(x);% 运行MCMC[theta_chain, logp](@ref)= run_mcmc(target, theta_init, num_iter);% 后处理theta_samples = theta_chain(burn_in:end,:);plot_results(theta_samples);
endfunction [chain, logp](@ref)= run_mcmc(target, theta0, num_iter)% 初始化d = length(theta0);chain = zeros(d,num_iter);chain(:,1) = theta0;logp = zeros(num_iter,1);% 初始对数概率logp(1) = target(theta0);% 迭代采样for i = 2:num_iter% 生成候选点theta_star = propose_step(chain(:,i-1));% 计算对数概率差logp_star = target(theta_star);log_alpha = logp_star - logp(i-1);% 接受/拒绝if log(rand) < log_alphachain(:,i) = theta_star;logp(i) = logp_star;elsechain(:,i) = chain(:,i-1);logp(i) = logp(i-1);endend
end
五、性能优化
- 预处理数据：对高维数据进行降维处理（PCA/T-SNE）
- 稀疏化采样：在低概率区域减少采样密度
- GPU加速：使用
gpuArray加速矩阵运算 - 诊断监控：实时绘制接受率与能量曲线
六、资源推荐
- 经典文献： Geman S., Geman D. (1984) 《Stochastic Relaxation》 Neal R. M. (2011) 《Handbook of Markov Chain Monte Carlo》
- 代码：多能源系统优化包含粒子群SVM等优化算法 www.youwenfan.com/contentcnm/80721.html
- MATLAB工具箱： Statistics and Machine Learning Toolbox Global Optimization Toolbox