简介基于MCMC马尔科夫-蒙特卡洛抽样的MATLAB仿真资源面向本硕博及教研人员适用于学习Metropolis-Hastings抽样、贝叶斯推断与概率分布采样等典型应用场景。资源包共5个文件包含3个MATLAB脚本、1个说明文档和1个avi操作录像压缩包大小仅753KB轻量易部署。脚本由主程序与两个子函数构成主程序与子函数分离便于逐段理解马尔科夫链状态转移、接受拒绝准则与积分估计的实现流程配套txt说明与操作录像则演示了在MATLAB2021a及以上版本中的正确启动方式——需将当前文件夹切换至工程路径并运行主文件而非直接运行子函数可帮助使用者避开常见运行错误。目前已有2381人学习下载适合作为MCMC算法课程设计、论文仿真或自学入门的参考资料借助实操视频可快速复现抽样过程并验证算法效果。1. 基于MCMC马尔科夫-蒙特卡洛抽样matlab仿真到底在解什么问题做贝叶斯推断的人迟早会撞上一堵墙后验分布写得出公式却算不出数值。抛硬币的二项分布还能手推共轭先验一旦模型变成层次结构、混合分布或者带隐变量的状态空间后验的归一化常数就是一个高维积分解析解不存在数值积分在维数超过四五个时直接失效。MCMC就是把“算积分”换成“抽样本”——让马尔可夫链的稳态分布恰好等于目标后验链走出来的样本就是对后验的近似描述。用MATLAB做这件事优势在于矩阵化计算、内置概率分布对象、以及调试可视化的一体化环境。这篇博客的落地路径是先讲清抽样原理和链的收敛条件再给出在MATLAB里可复现的Metropolis-Hastings和Gibbs采样实现接着解决“仿真发散”和“链不收敛”的排查方法最后用MATLAB把采样过程导出成操作视频方便汇报和教学。2. 马尔科夫-蒙特卡洛的核心蒙特卡洛积分、细致平衡与转移核2.1 蒙特卡洛采样为什么能替代解析积分蒙特卡洛的核心思想是用样本均值替代积分期望。假设想计算函数g(x)关于分布p(x)的期望解析做法是求∫g(x)p(x)dx但高维情况下这个积分无法计算。蒙特卡洛的做法是从p(x)中抽取N个独立样本用(1/N)Σg(x_i)近似期望。大数定律保证N足够大时近似收敛中心极限定理则给出了误差的收敛速度——大概是N的平方根分之一。这意味着精度每提高一个数量级样本量要增加一百倍所以单纯增加抽样次数不是好办法真正要做的是让每一个样本都来自目标分布。MCMC与普通蒙特卡洛的区别在于样本并非独立。它构造一条马尔可夫链当前状态只依赖上一个状态通过精心设计的转移概率让链在足够长的时间后收敛到目标分布。链在平稳状态下的样本虽然相关但依然可以用来做期望估计。这就是为什么MCMC适用于高维、非标准、无法直接采样的分布而普通蒙特卡洛只能用于已知简单分布。2.2 细致平衡条件与转移核的设计约束要让一条马尔可夫链收敛到目标分布p(x)必须满足细致平衡条件对于任意两个状态x和x有p(x)T(x→x) p(x)T(x→x)。T表示转移核即从当前状态跳转到另一个状态的概率密度。满足细致平衡能保证p(x)是链的不变分布但不保证链一定收敛还要求链具有不可约性和非周期性。不可约指从任意状态都能到达其他状态非周期指不会陷入固定步长的循环振荡。用MATLAB理解这些概念最直观的方式是先不看MCMC而是手动构造一个离散状态的马尔可夫链。% 定义一个3状态马尔可夫链的转移矩阵 P [0.7 0.2 0.1; 0.3 0.4 0.3; 0.1 0.3 0.6]; % 初始状态分布从状态1开始 pi0 [1 0 0]; pi_current pi0; chain_states zeros(2000, 1); for step 1:2000 % 按当前状态所在行进行随机跳转 chain_states(step) find(mnrnd(1, P(pi_current 1, :))); pi_current zeros(1, 3); pi_current(chain_states(step)) 1; end % 计算各状态出现频率对比理论平稳分布 empirical_pi histcounts(chain_states, 1:4) / 2000; [eigen_vec, ~] eig(P); stationary eigen_vec(:,1) / sum(eigen_vec(:,1)); disp(经验频率:); disp(empirical_pi); disp(理论平稳分布:); disp(real(stationary));这段代码从状态1出发按转移矩阵P反复跳转2000步。mnrnd是MATLAB内置的多项分布随机数生成函数第一个参数1表示每次抽一个样本第二个参数是当前状态对应的转移概率行。迭代结束后统计每个状态出现的频次并与P矩阵的特征向量对比——平稳分布就是特征值为1对应的左特征向量归一化结果。这个最小示例用来验证直观感受不管初值如何链跑足够长后状态访问频率总会稳定在固定比例上。MCMC的原理与此一样区别只在转移矩阵是连续状态下的核函数且目标分布不是固定矩阵的特征向量而是待抽样的后验。2.3 随机游走转移核与接受-拒绝框架连续状态MCMC中最常见的构造方式是随机游走Metropolis-Hastings算法简称MH。其做法是当前状态记作θ从提议分布q(θ|θ)中生成候选状态θ然后计算接受概率α min(1, p(θ)q(θ|θ) / (p(θ)q(θ|θ)))以概率α接受候选状态作为链的下一个状态否则保持原状态。提议分布的选择直接影响采样效率。高斯分布是应用最广的提议分布θ θ σ·ε其中ε服从标准正态分布。对称提议分布下q(θ|θ) q(θ|θ)接受率退化为后验密度比值。σ步长太小接受率高但链移动缓慢样本强相关有效样本量低σ步长太大候选点经常落到后验的尾巴区域被拒绝的频率变高链长时间停留原地。经验法则是调节σ让接受率落在0.2到0.5之间高斯目标变量维度较低时0.44附近比较合适。3. 在MATLAB中实现Metropolis-Hastings抽样目标分布、候选分布与参数整定3.1 用MH采样双峰高斯混合分布的完整代码真实建模中后验经常是多峰的这里用一个双峰高斯混合作为目标分布展示MCMC如何从多峰分布中正确抽样。% 目标分布的双峰高斯混合 mu [-3; 3]; sigma [1.0; 1.2]; weight [0.4; 0.6]; log_target (x) log(weight(1)*normpdf(x, mu(1), sigma(1)) ... weight(2)*normpdf(x, mu(2), sigma(2))); % MH参数 n_samples 20000; % 总采样数 burn_in 5000; % 燃烧期样本数 sigma_prop 2.0; % 提议分布标准差 theta0 -5; % 初始值 theta zeros(n_samples, 1); theta(1) theta0; n_accept 0; % 随机游走MH主循环 for t 2:n_samples % 从高斯提议分布生成候选点 theta_cand theta(t-1) sigma_prop * randn(); % 计算对数接受率用log避免数值下溢 log_alpha log_target(theta_cand) - log_target(theta(t-1)); log_alpha min(0, log_alpha); % log(min(1, ratio)) % 接受-拒绝判定 if log(rand()) log_alpha theta(t) theta_cand; n_accept n_accept 1; else theta(t) theta(t-1); end end % 输出接受率并绘图 accept_rate n_accept / (n_samples - 1); fprintf(接受率: %.3f\n, accept_rate); samples theta(burn_in 1:end); figure; histogram(samples, 200, Normalization, pdf); hold on; x_range linspace(-8, 8, 500); p_true weight(1)*normpdf(x_range, mu(1), sigma(1)) ... weight(2)*normpdf(x_range, mu(2), sigma(2)); plot(x_range, p_true, r-, LineWidth, 2); legend(MCMC样本直方图, 真实目标分布);代码里的关键点是全程使用对数密度。目标分布是混合高斯时密度值本身可能极小直接计算比值会长期等于0或inf采样被迫中断。对数化之后比值变为减法log(rand())与log_alpha的比较在数值上稳定得多。建议参数含义如下sigma_prop是提议分布标准差决定候选状态的平均跳跃距离burn_in是燃烧期样本这段样本尚未收敛到稳态直接丢弃不参与统计n_accept用于输出真实接受率与预设经验值对照。运行后两个峰值都会被采出来但链条在两个峰之间切换的频率取决于sigma_prop——步长太小会长期困在单一峰内。3.2 对数概率比的计算为什么比直接除法可靠直接计算接受率α min(1, p(θ)/p(θ))在概念上简单实际执行时却受限于浮点数精度。后验密度通常涉及连乘的似然项比如观测数据有1000个点每个点的似然在0.1到0.9之间连乘之后数值直接下溢为0。0除以0没有意义。取对数之后连乘变成求和密度值的大小范围从10的负几百次方缩放到几百以内完全落入双精度浮点的表示区间。% 直接计算与对数计算的对比示例 x 1:1000; likelihood_vec 0.5 * exp(-0.5 * (x - 500).^2 / 100^2); p_direct prod(likelihood_vec); % 下溢为0 p_log sum(log(likelihood_vec)); % 对数域正常 fprintf(直接乘积: %g\n, p_direct); fprintf(对数求和: %.2f\n, p_log);这段代码模拟了1000个观测点的似然连乘。prod运算结果在MATLAB中显示为0因为中间步骤已经低于最小浮点数。而sum(log(...))得到一个负几百的数值完全可用。MH算法里只需要后验密度的比值对数空间下log(p(θ))与log(p(θ))相减即可得到log_alpha不需要恢复原始尺度。如果模型包含先验分布把先验的对数密度加上似然的对数密度整个计算过程始终保持在对数域。3.3 提议步长sigma的整定策略与接受率对照表sigma_prop取值接受率区间链的行为表现适用场景0.1 ~ 0.30.7~0.9几乎每次接受但相邻样本高度相关有效样本量极低仅适合低维且目标分布方差很小的情况0.5 ~ 1.50.4~0.6移动效率较高自相关适中最常用区间单峰且尺度接近标准正态的后验2.0 ~ 4.00.2~0.35频繁拒绝但跳跃幅度大链混合速度尚可双峰分布或多峰时适应能力较强5.0 以上低于0.1链长时间停在原处极少数大跳跃导致严重偏态不推荐出现时需缩小步长接受率并不是越高越好。接受率接近1说明提议分布的标准差远小于目标分布的标准差链在局部范围内转圈遍历整个后验需要极长的链长。接受率过低说明提议的跳跃经常落在后验支撑集外面浪费了计算资源。一个实际做法是先用短链做预实验以1000步为一次迭代调sigma观察接受率落在0.2~0.5区间后再正式跑长链。对高维问题单一随机游走步长会失效此时应当换用自适应MH算法或者按维度分组采样。4. 仿真发散与链不收敛MATLAB下的诊断方法和处理策略4.1 仿真发散的具体表现轨迹图、直方图与自相关仿真发散在MCMC语境里指链没有收敛到目标后验或者在两个模态之间切换不足统计结果失真。肉眼可见的标志是轨迹图出现长平台段——链长期停在同一数值附近不动然后突然跳到另一片区域。另一个标志是直方图与目标分布明显偏离比如双峰分布只采到一个峰。自相关图衰减缓慢则说明相邻样本之间信息冗余严重有效样本量不足。在MATLAB里写一个简单的诊断脚本可以让这些问题当场暴露。samples_good theta(burn_in 1:end); % 轨迹图观察链是否遍历整个支撑区间 figure; subplot(2,1,1); plot(samples_good, LineWidth, 0.3); title(轨迹图); xlabel(迭代次数); ylabel(采样值); % 自相关图观察链的记忆长度 subplot(2,1,2); [autocorr_lags, autocorr_vals] autocorr(samples_good, NumLags, 100); stem(autocorr_lags, autocorr_vals, MarkerSize, 3); title(sprintf(自相关函数滞后期100)); xlabel(lag); ylabel(自相关系数);轨迹图如果像毛毛虫那样围绕一条水平线均匀抖动说明链平稳。如果出现明显断层或者长时间上升/下降趋势说明链还处于燃烧期需要丢弃更多前期样本或增大burn_in。自相关函数应该在若干滞后之后迅速跌向0如果滞后50之后仍然高于0.2说明相邻样本高度相关应该考虑增大sigma_prop或对样本做稀疏化存储每隔几步记录一个样本。4.2 39个潜在尺度缩减因子R_hat用多链判断是否收敛单链轨迹图可能产生误判——链在一个局部模态停留足够久也会呈现“看似平稳”的画面。业界通行做法是同时跑多条链比较链间方差与链内方差这就是Gelman-Rubin诊断也叫R_hat。收敛标准是R_hat小于1.1越接近1越好。MATLAB中虽然没有直接封装R_hat计算但用统计工具箱可以快速手写。n_chains 4; n_iter 15000; burn_in 5000; chain_samples zeros(n_iter, n_chains); % 并行跑4条链初始值分散在不同位置 start_points [-6; -1; 2; 6]; parfor c 1:n_chains chain_samples(:, c) run_mh_sampler(start_points(c), n_iter); end % 计算R_hat samples_post chain_samples(burn_in 1:end, :); n size(samples_post, 1); chain_means mean(samples_post, 1); grand_mean mean(chain_means, 2); B n / (n_chains - 1) * sum((chain_means - grand_mean).^2); W mean(var(samples_post, 0, 1), 2); var_est (n - 1) / n * W B / n; R_hat sqrt(var_est / W); fprintf(R_hat %.3f\n, R_hat);multiple chains的核心意义在于暴露局部模态问题。如果所有链初始值分散在参数空间不同区域最终都收敛到几乎相同的均值与方差那可以确信链真正遍历了后验。R_hat大于1.1时不能信任当前采样结果常见解决办法是延长燃烧期、增大链总数或改用自适应提议分布。4.3 处理多峰分布的收敛陷阱模拟退火与并行回火多峰后验是MCMC仿真发散的重灾区。普通MH从某个峰附近出发后要跨越两个峰之间的低概率区域才能到达另一个峰而低概率区域的接受概率极低链几乎不可能完成跨越。解决办法分两类一类是把目标分布取幂p_α(θ) ∝ p(θ)^(1/T)温度T越高分布越平坦峰谷之间的差异被抹平链可以自由巡游整个空间随后逐渐降低T恢复到原始分布。这就是模拟退火。另一类是并行回火同时维护多条不同温度的链每隔一定步数按Metropolis准则交换两条链的状态高温链帮低温链跳出局部陷阱。MATLAB实现并行回火的关键代码片段function swapped_samples parallel_temper_mh(log_target, n_chains, n_iter, swap_interval, temps) n_samples n_iter * n_chains; chain_states randn(n_chains, 1) * 5; swapped_samples zeros(n_samples, 1); idx 1; for iter 1:n_iter for c 1:n_chains proposal chain_states(c) 1.5 * randn(); log_alpha (log_target(proposal) - log_target(chain_states(c))) / temps(c); if log(rand()) min(0, log_alpha) chain_states(c) proposal; end end % 相邻链状态交换 if mod(iter, swap_interval) 0 for c 1:(n_chains - 1) log_swap (log_target(chain_states(c 1)) - log_target(chain_states(c))) * ... (1 / temps(c) - 1 / temps(c 1)); if log(rand()) min(0, log_swap) temp_state chain_states(c); chain_states(c) chain_states(c 1); chain_states(c 1) temp_state; end end end for c 1:n_chains swapped_samples(idx) chain_states(c); idx idx 1; end end end温度序列一般按几何级数排布最低温为1最高温让链能自由跨越两个峰。链间交换间隔不宜太频繁否则状态交换开销过大也不宜太稀疏否则高温链无法帮助低温链。实际运行时先试swap_interval 50观察低温链的轨迹图是否出现模态切换再逐步调整。5. 把MCMC采样过程做成操作视频MATLAB动画录制与结果导出5.1 用animatedline动态绘制轨迹和高斯核密度估计汇报和教学场景里静态直方图不足以展示MCMC的抽样过程。动态呈现采样点一波一波落入直方图观众才能直观理解马尔可夫链逐步逼近目标分布的过程。MATLAB的animatedline函数可以在循环中不断追加数据点并实时刷新图形搭配drawnow即可形成流畅动画。% 动态绘制MCMC采样过程动画 figure; h_hist histogram(BinEdges, -8:0.2:8, BinCounts, zeros(1, 80), ... FaceColor, [0.3 0.6 0.9], EdgeColor, none); hold on; % 目标分布参考曲线 x_line linspace(-8, 8, 300); p_line weight(1)*normpdf(x_line, mu(1), sigma(1)) ... weight(2)*normpdf(x_line, mu(2), sigma(2)); plot(x_line, p_line, r-, LineWidth, 2); ylim([0 0.3]); xlabel(θ); ylabel(概率密度); % 动态更新直方图 bin_edges -8:0.2:8; hist_counts zeros(1, length(bin_edges) - 1); total_draws 0; for t 1:2000 % theta_precomputed为事先跑好的MH链 idx find(theta_precomputed(t) bin_edges(1:end-1) ... theta_precomputed(t) bin_edges(2:end), 1); if ~isempty(idx) hist_counts(idx) hist_counts(idx) 1; end total_draws total_draws 1; if mod(t, 5) 0 set(h_hist, BinCounts, hist_counts / total_draws / 0.2); title(sprintf(MCMC采样过程 | 第 %d 次迭代, t)); drawnow; end end这个动画展示的是归一化直方图逐步逼近目标分布曲线的过程。每5步刷新一次可以降低渲染开销2000次迭代的动画大约5秒内完成。运行前需要预先用MH算法生成theta_precomputed序列动画只是将采样过程重放一遍避免动画播放期间的回放与真实采样相互干扰。5.2 输出为GIF或MP4录制操作视频与分享格式选择MATLAB导出动画有两种推荐方式。一是输出为GIF适合插入PPT或网页展示二是输出为MP4适合完整操作视频录制。转换逻辑都是先把动画每一帧捕获下来再编码成目标格式。% 方案一导出GIF figure; gif_filename mcmc_sampling.gif; for t 1:500 % 更新图形的代码与5.1节相同 % ... frame getframe(gcf); [A, map] rgb2ind(frame.cdata, 256); if t 1 imwrite(A, map, gif_filename, gif, LoopCount, Inf, DelayTime, 0.05); else imwrite(A, map, gif_filename, gif, WriteMode, append, DelayTime, 0.05); end end % 方案二导出MP4视频 video_file VideoWriter(mcmc_sampling.mp4, MPEG-4); video_file.FrameRate 20; open(video_file); for t 1:500 % 更新图形的代码 % ... frame getframe(gcf); writeVideo(video_file, frame); end close(video_file);GIF的优势是到处能播放缺点是有256色限制直方图色彩过渡会出现轻微色带。MP4色彩保真度更高文件体积也更小但需要播放器支持。DelayTime参数控制每帧停留时间0.05秒对应20帧每秒的输出速度符合人眼视觉暂留特性。导出视频之前建议把figure窗口置于屏幕最前并固定尺寸getframe捕获的画面大小与figure窗口显示尺寸严格一致窗口大小影响视频分辨率。5.3 操作视频中的关键标注高斯核密度叠加与接受率显示单纯动画还不够操作视频里最好把每次迭代是“接受”还是“拒绝”用颜色标出来接受的点用深色落进直方图拒绝时当前帧高亮显示并保留原状态。另一种有效标注是动态显示累计接受率让观看者看到链的调整过程。% 在绘图区域左上角动态显示接受率文本 text_handle annotation(textbox, [0.15 0.8 0.3 0.1], ... String, , EdgeColor, none, ... FontSize, 12, FontWeight, bold); for t nanot~:2000 % ...绘图更新代码... current_accept_rate n_accept / t; set(text_handle, String, sprintf(累计接受率: %.2f%%, current_accept_rate * 100)); drawnow; end这段代码用annotation在图形坐标系的相对位置绘制文本不受坐标轴范围变化影响。接受率文本随着迭代进行逐渐稳定到目标区间观看者能看到参数整定的动态过程。最后补充一句做操作视频时把燃烧期用红色竖线在轨迹图上标出来观众一眼就能分辨哪些样本进入统计、哪些被丢弃。本文还有配套的精品资源点击获取