简介基于MCMC马尔科夫-蒙特卡洛抽样的Matlab仿真资源面向本硕博及科研教学人群适合用于马尔科夫链蒙特卡洛抽样算法的编程学习、教学演示与算法验证。资源共5个文件以Matlab源程序为主3个.m另含1个avi操作录像和1个txt说明文档压缩包整体仅753KB轻量紧凑、即下即用。代码层面包含主运行脚本、被调用的子函数以及MCMC抽样算法实现可直接在Matlab 2021a及以上版本运行操作录像展示从路径设置、脚本执行到结果观察的完整过程txt说明则补充了主程序与子函数的运行关系及注意事项学习者可按视频快速复现抽样过程。目前已积累2381人学习下载适合希望通过源码加演示系统理解马尔科夫链、建议分布与接受-拒绝机制的初学者也便于教师在课堂上作为实验案例使用。1. 你以为MCMC是黑箱其实它就是一条能收敛的马尔科夫链当后验分布的归一化常数求不出来、解析解不存在时贝叶斯推断就卡住了。MCMC把问题从“求积分”换成“抽样”构造一条以目标分布为平稳分布的马尔科夫链沿着链子走下去轨迹的样本分布自然逼近那个求不出来的后验期望值用样本均值替代。这条思路看着绕工程上却非常可靠。matlab仿真里最常用的两个入口是Metropolis-HastingsMH和Gibbs采样器前者只需要未归一化的密度表达后者把联合后验拆成一组好采的条件分布。这篇内容面向需要用MATLAB把MCMC跑通、又不想只按工具箱黑箱调参的人从平稳分布原理讲到接受率、提议方差、收敛诊断和贝叶斯回归实战代码可直接复制运行。2. MCMC的马尔科夫链机制平稳分布、细致平衡与混合效率2.1 从“求不出积分”到“抽得出样本”蒙特卡洛在数据增强里的角色后验推断的通用形式是算后验期望 E[f(θ)|y] ∫f(θ)π(θ|y)dθ其中π(θ|y) ∝ L(y|θ)π(θ)。分子好写分母是那个归一化常数∫L(y|θ)π(θ)dθ维度一高就积不出来。蒙特卡洛的想法是从π(θ|y)直接抽取大量样本均值收敛到上述期望。问题在于“直接抽”在很多分布上做不到——逆变换要求可积的CDF拒绝采样在高维空间接受率指数级下降重要性采样在重尾时方差爆炸。MCMC换了个思路构造一条转移链让它在长期运行中“生成”目标分布的样本不需要求归一化常数也不需要直接采样。这条链的关键不在样本是否独立而在转移核P(x→x′)能使目标分布π在每一步之后都保持不变也就是不动点条件π πP。工程上更常用的是一个充分条件叫detailed balance细致平衡π(x)P(x→x′) π(x′)P(x′→x)满足它π一定是该链的平稳分布。MH算法的全部工作就是设计转移核让这个条件成立给定当前x先从建议分布q(x′|x)抽候选再以接受率决定是否移动。2.2 转移核与平稳分布MH如何满足detailed balanceMH的接受率形式是 α(x, x′) min(1, π(x′)q(x|x′) / (π(x)q(x′|x)))。如果提议分布对称比如高斯随机游走提议x′ x N(0, σ²)q(x|x′) q(x′|x)可以约掉接受率退化为min(1, π(x′)/π(x))这就是Metropolis采样。这个简化是matlab仿真里最常用的结构写目标分布时只需要定义联合密度而不用管归一化常数因为比值里常数消掉了。这正是MCMC最大的工程红利——很多应用问题里写得出似然×先验但归一化常数根本没闭式。接受率高低成了转移核设计里的核心变量。接受率过高意味着新样本离当前点很近链在局部徘徊混合慢过低意味着每步都在拒绝链原地踏步。高维情况下理论最优接受率在0.234附近Roberts等人对高斯目标的渐近结果一维目标大约0.44。仿真中不必精确逼近这个数只需监控接受率落在0.2~0.5区间偏离太多就调整提议方差。混合效率直接决定后面看到的自相关强度。2.2.1 接受率公式的对称性简化与数值风险细致平衡是局部条件它确保任意一步转移中从x到x′的期望流等于反向流。把两侧对x积分会发现π在转移作用下不变。这个证明在MATLAB中不必显式写出但它决定了接受率的结构分子分母必须包含前向概率π(x′)q(x′|x)和反向概率π(x)q(x|x′)二者比值失衡时通过随机接受修正。很多仿真发散的问题追根溯源是写接受率时把q项当成常量约掉破坏了详细平衡条件。2.3 burn-in、混合与自相关样本不是独立同分布的从任意初始值出发链需要一段“预热”才能进入平稳区域这一段称为burn-in样本要丢弃。进入平稳后相邻样本高度相关——这是马尔科夫链的固有属性不是bug。相关越强链的有效样本量越小。量化上使用自相关长度τ有效样本量ESS N/ττ越大信息越少。仿真时用trace plot横轴迭代次数纵轴样本值直接看混合状态理想情况下链在两个模式之间频繁穿梭而不是长时间困在一侧。三张图画完MCMC仿真才算真的“可见”了。3. 用MATLAB从零写Metropolis-Hastings采样器对准双峰混合高斯3.1 目标分布与提议分布的选型为什么用随机游走高斯提议双峰混合高斯是MCMC最常用的试验场分布有两个峰值如果链只在一个峰附近徘徊trace plot一眼就能看出来。设定目标分布为π(x) 0.6·N(x; -2, 0.8²) 0.4·N(x; 3, 1.2²)。这个分布没有解析采样方法但密度表达式简单。提议分布采用对称高斯随机游走x′ x σ·randn()。对称性的好处是接受率里没有q项实现简单缺点是高维时固定方差难以匹配不同维度尺度后面自适应MCMC章节会解决。在确定性仿真场景下σ的取值就是最重要的调参旋钮。先定义对数密度函数。这里有一个资深用户常踩的坑不要直接用exp相加再取log两个高斯分量强度差异大时exp(w)会下溢为0log变成-Inf。应该在log域计算用log-sum-exp技巧function logp log_target(x) % 双峰高斯混合的对数密度返回标量 % 每个分量的log权重 log高斯核归一化常数在MH比值中抵消 logw zeros(size(x, 1), 2); % 第一个分量均值-2标准差0.8权重0.6 logw(:, 1) log(0.6) - 0.5 * ((x 2) / 0.8).^2; % 第二个分量均值3标准差1.2权重0.4 logw(:, 2) log(0.4) - 0.5 * ((x - 3) / 1.2).^2; % log-sum-exp先减最大值再取exp防止指数下溢 m max(logw, [], 2); logp m log(sum(exp(logw - m), 2)); end这段函数的关键在于log-sum-exp把公共最大值m提出来exp(logw - m)的值域被控制在(0,1]log(sum)不会爆炸。均值、标准差、权重是目标分布的物理意义参数而- log(0.8·sqrt(2π))这类常数可以放心丢掉因为MH接受率是比值。3.2 主循环代码与接受率log域计算避免下溢MH主循环如下。以x00为初始状态前5000次作为burn-in丢弃后20000次保留% 主循环参数 N 20000; % 采样总数 burnin 5000; % 预热步数 sigma 1.5; % 提议标准差随机游走步长 x 0; % 初始状态 samples zeros(N, 1); acc 0; % 接受计数器 for t 1:(N burnin) % 高斯随机游走提议对称分布q项在MH中抵消 x_star x sigma * randn(); % log域接受率接受概率 exp(logp(x_star) - logp(x)) log_alpha log_target(x_star) - log_target(x); if log(rand()) log_alpha x x_star; acc acc 1; end % 丢弃burn-in阶段其余样本保存 if t burnin samples(t - burnin) x; end end acceptance_rate acc / (N burnin); fprintf(接受率: %.3f\n, acceptance_rate);主线是“提议-接受/拒绝-记录”。log(rand()) log_alpha把接受概率转化为一次均匀随机数与阈值的比较避免了对alpha再取exp的数值风险。计数器放在if内部分母包含burnin因为预热阶段的接受率同样反映提议方差是否合适。运行结束后接受率约在0.3附近说明步长合理低于0.1要减小sigma高于0.6则增大sigma。3.2.1 提议方差σ怎么设三个档位与接受率对照sigma的选择直接决定采样质量经验上分成三个档位sigma取值接受率区间自相关表现典型问题0.1~0.30.7相邻样本高度相关ESS小链移动太慢burn-in极长0.8~2.00.2~0.5适中trace快速穿梭最理想区间5以上0.05链长时间停在原地大量拒绝仿真发散感强提示自相关强不是采样失败而是有效样本量不足的信号。此时优先调sigma而不是盲目增加总迭代数。3.3 采样结果的三个验证图trace、直方图、自相关画验证图的代码集中在一次figure里subplot(2, 2, 1); plot(1:N, samples); title(Trace Plot); xlabel(迭代次数); ylabel(x); subplot(2, 2, 2); histogram(samples, 60, Normalization, pdf); hold on; % 把真实混合分布画上去做对比 xs linspace(-6, 7, 400); true_pdf 0.6 * normpdf(xs, -2, 0.8) 0.4 * normpdf(xs, 3, 1.2); plot(xs, true_pdf, r-, LineWidth, 2); legend(MCMC样本, 真实密度); subplot(2, 2, 3); autocorr(samples, 50); % 自相关图横轴lag纵轴自相关 fprintf(样本均值: %.3f\n, mean(samples)); fprintf(样本方差: %.3f\n, var(samples));3.3.1 三个图分别回答什么问题trace plot回答“链收敛了吗”。收敛的链纵向抖动幅度稳定两个峰之间周期性跨越若长时间只在一个峰附近要么是burn-in不够要么是sigma太小。直方图回答“样本分布对吗”。MCMC样本的直方图与真实密度红色曲线贴合说明抽样系统正确贴合但形状粗糙说明样本数不够或自相关太高。自相关图回答“信息量真的够吗”。自相关从lag0的1.0衰减到0的速度越快越好若滞后30仍然高于0.5有效样本量会显著低于名义样本量。操作视频的演示里通常把这三张图并列放在同一个figure窗口运行结束一眼确认收敛状态。下一步顺着变异方向走当变量维度增多靠手调sigma不现实就要引入Gibbs采样和自动协方差调节。4. 从Gibbs到自适应MCMC样本效率、R-hat与Geweke诊断4.1 Gibbs采样用满条件分布更新每一个分量当目标分布是多维联合分布而每个分量的满条件分布有闭型时Gibbs采样比MH更高效——不需要调节提议方差也没有拒绝步。做法是把θ拆成(θ1, …, θd)每一轮依次从p(θi|θ_{-i}, y)采样替换当前值。这个循环保持每个条件分布为平稳分布组合起来就采样了联合后验。工程上最常见的场景是共轭先验下的贝叶斯模型正态-正态、正态-逆伽马组合都有闭式条件分布。MATLAB中实现Gibbs的骨架如下% Gibbs采样骨架三变量目标条件分布假设均有闭型 theta zeros(3, 1); n_iter 10000; samples_gibbs zeros(n_iter, 3); % 伪代码示意实际条件分布需按模型代入 for iter 1:n_iter % 依次更新每个分量直接以更新的值作为后续分量条件 theta(1) conditional_sample_1(theta(2), theta(3)); theta(2) conditional_sample_2(theta(1), theta(3)); theta(3) conditional_sample_3(theta(1), theta(2)); samples_gibbs(iter, :) theta; end更新顺序有讲究第1个分量总是使用最新的θ2、θ3第2个分量使用刚更新的θ1和旧的θ3。这种“始终用新值”的扫描方式称为systematic scan是MATLAB中条件更新矩阵运算的标准写法能加快混合速度。4.2 自适应提议2.38²/d缩放与协方差学习期Gibbs要求条件分布可采样但很多后验并不满足。回到MH框架在多维情况下把单变量随机游走升级为多元高斯提议x′ x N(0, σ²·Σ)其中Σ是目标分布的协方差估计σ是缩放因子。常用的经验法则是σ 2.38 / sqrt(d)d为目标维度这个缩放因子在高斯目标下渐近最优。MATLAB里通过chol分解生成多元高斯样本% 自适应提议的协方差学习过程 d 2; Sigma_adapt eye(d) * 0.02; % 初始小协方差 L_adapt chol(Sigma_adapt, lower); N_total 8000; learn_period 3000; % 学习期之后提议固定 x zeros(d, 1); samples_adapt zeros(N_total, d); for t 1:N_total % 用当前协方差矩阵的Cholesky因子生成多元高斯增量 increment L_adapt * randn(d, 1); x_star x 2.38 / sqrt(d) * increment; log_alpha log_target_vec(x_star) - log_target_vec(x); if log(rand()) log_alpha x x_star; end samples_adapt(t, :) x; % 仅在学习期内更新协方差估计学习期结束后固定提议 if t learn_period mod(t, 50) 0 Sigma_adapt cov(samples_adapt(max(1, t-500):t, :)); Sigma_adapt Sigma_adapt 1e-8 * eye(d); % 正则化防奇异 L_adapt chol(Sigma_adapt, lower); end end这里两个必调参数学习期learn_period和正则项系数1e-8。学习期必须小于总迭代数结束后提议协方差固定否则链的非马尔科夫性破坏收敛理论。正则化防止协方差矩阵在样本少时奇异导致chol失败。注意cov()接收的是过去500个样本形成的窗口滑动窗口比全历史更能捕捉局部协方差结构。提示自适应MCMC只是在学习期用过去的样本来“调参”正式采样仍要求提议不变。仿真中若看到学习期之后接受率骤降多半是learn_period设置过短协方差还没收敛就锁定了。4.3 收敛诊断三件套R-hat、有效样本量与Geweke检验把MCMC当黑箱用时收敛诊断是保险。R-hat需要至少两条独立链比较链间方差与链内方差function Rhat compute_rhat(chains) % chains: M条链 × N次迭代矩阵 [M, N] size(chains); % 链均值的方差 → 链间方差B chain_means mean(chains, 2); overall_mean mean(chain_means); B N / (M - 1) * sum((chain_means - overall_mean).^2); % 每条链内方差的平均 → 链内方差W W mean(var(chains, 0, 2)); % 总方差估计混合公式 var_est (N - 1) / N * W B / N; Rhat sqrt(var_est / W); end计算时注意var(x, 0, 2)的维度参数0表示除以N-12表示按行计算。R-hat越接近1越好。有效样本量ESS针对单条链的自相关结构function ess effective_sample_size(x) % 通过自相关估计自相关长度tau再得到ESS n numel(x); xc xcorr(x - mean(x), coeff); rho xc(n:end); % 保留非负lag一侧 % 截断自相关尾部噪声对长lag求和影响大 max_lag min(1000, n - 1); tau 1 2 * sum(rho(2:max_lag 1)); ess n / tau; end自相关尾部噪声会让tau被高估max_lag上限正是为了压制这一点。Geweke检验取前10%与后50%的样本均值比较构造的z统计量落在[-1.96, 1.96]内表示未检测到收敛异常。三者的用法对照指标阈值/参考值结果解读R-hat1.1多链一致可认为已收敛ESS400后验均值估计稳定Geweke z|z| 1.96前后段均值无显著差异仿真发散时先看这三项再决定是加长学习期、增大迭代数还是回头检查接受率和提议方差。5. 实战闭环用MCMC拟合贝叶斯线性回归并输出后验可信区间5.1 模型设定与条件后验推导把前面的技术落到最常见的回归场景。模型为y Xβ εε ∼ N(0, σ²I)。配共轭先验β ∼ N(β0, Λ0⁻¹)σ² ∼ IG(a0, b0)。该设定下条件后验有闭型β|σ², y ∼ N( (XᵀX/σ² Λ0)⁻¹(Xᵀy/σ² Λ0β0), (XᵀX/σ² Λ0)⁻¹ )σ²|β, y ∼ IG( a0 n/2, b0 (y - Xβ)ᵀ(y - Xβ)/2 )。第二个条件分布只依赖残差这是Gibbs能闭环的关键。生成仿真数据rng(42); n 200; X [ones(n, 1), randn(n, 1), randn(n, 1)]; % 截距 两个自变量 beta_true [0.5; -1.2; 0.8]; y X * beta_true 0.5 * randn(n, 1);5.2 Gibbs循环的MATLAB实现按条件更新β和σ²% 共轭先验参数 d size(X, 2); beta0 zeros(d, 1); Lambda0 eye(d) * 0.01; a0 0.01; b0 0.01; % 初始值 beta beta0; sigma2 1.0; n_iter 12000; burnin 2000; beta_samples zeros(n_iter - burnin, d); sigma2_samples zeros(n_iter - burnin, 1); for iter 1:n_iter % 1. 更新beta多元正态协方差矩阵来自XX/σ² Λ的逆 Xbar X * X / sigma2 Lambda0; beta_mean Xbar \ (X * y / sigma2 Lambda0 * beta0); Sigma_beta inv(Xbar); Lb chol(Sigma_beta, lower); beta beta_mean Lb * randn(d, 1); % 2. 更新sigma²逆伽马注意gamrnd第二个参数是scale不是rate resid y - X * beta; a_post a0 n / 2; b_post b0 0.5 * (resid * resid); sigma2 1 / gamrnd(a_post, 1 / b_post); if iter burnin beta_samples(iter - burnin, :) beta; sigma2_samples(iter - burnin) sigma2; end end两个常见的坑都在这里gamrnd(a, b)第二个参数是scale逆伽马的rate参数需要取倒数所以写1/b_postbeta_mean用矩阵左除Xbar \ (...)而不是inv(Xbar)*...数值稳定性更好尤其当Xbar接近病态时。5.3 输出后验结果并与OLS对比采样结束后直接计算后验统计量beta_mean_mcmc mean(beta_samples, 1); beta_ci quantile(beta_samples, [0.025, 0.975], 1); beta_ols (X * X) \ (X * y); fprintf(OLS: %.3f %.3f %.3f\n, beta_ols); fprintf(MCMC: %.3f %.3f %.3f\n, beta_mean_mcmc); fprintf(95%% CI:\n); disp(beta_ci);输出对照关系大致如下表实际运行会有小数位波动真实系数应落在CI内参数OLS估计MCMC后验均值95% CIβ0约0.47约0.48[0.28, 0.69]β1约-1.22约-1.22[-1.29, -1.15]β2约0.83约0.82[0.74, 0.89]这里的核心收益是“输出完整后验而不只是点估计”OLS给出点和标准误MCMC给出分布可以随时改置信水平、计算P(β1 0)这类问题。收尾时跑一下effective_sample_size(beta_samples(:,2))把ESS卡在400以上再汇报结果若低于阈值就增加迭代数而不是缩短burn-in。整套流程对照操作视频逐帧确认burn-in段恰好对应trace plot从初始值摆动到平稳区域的时刻协方差学习期结束对应接受率由震荡转平缓的转折点这两个视觉标志对齐了你手上这套基于MCMC马尔科夫-蒙特卡洛抽样matlab仿真才算真正跑通。本文还有配套的精品资源点击获取 用MATLAB搞定MCMC采样：MH算法、Gibbs与收敛诊断
简介基于MCMC马尔科夫-蒙特卡洛抽样的Matlab仿真资源面向本硕博及科研教学人群适合用于马尔科夫链蒙特卡洛抽样算法的编程学习、教学演示与算法验证。资源共5个文件以Matlab源程序为主3个.m另含1个avi操作录像和1个txt说明文档压缩包整体仅753KB轻量紧凑、即下即用。代码层面包含主运行脚本、被调用的子函数以及MCMC抽样算法实现可直接在Matlab 2021a及以上版本运行操作录像展示从路径设置、脚本执行到结果观察的完整过程txt说明则补充了主程序与子函数的运行关系及注意事项学习者可按视频快速复现抽样过程。目前已积累2381人学习下载适合希望通过源码加演示系统理解马尔科夫链、建议分布与接受-拒绝机制的初学者也便于教师在课堂上作为实验案例使用。1. 你以为MCMC是黑箱其实它就是一条能收敛的马尔科夫链当后验分布的归一化常数求不出来、解析解不存在时贝叶斯推断就卡住了。MCMC把问题从“求积分”换成“抽样”构造一条以目标分布为平稳分布的马尔科夫链沿着链子走下去轨迹的样本分布自然逼近那个求不出来的后验期望值用样本均值替代。这条思路看着绕工程上却非常可靠。matlab仿真里最常用的两个入口是Metropolis-HastingsMH和Gibbs采样器前者只需要未归一化的密度表达后者把联合后验拆成一组好采的条件分布。这篇内容面向需要用MATLAB把MCMC跑通、又不想只按工具箱黑箱调参的人从平稳分布原理讲到接受率、提议方差、收敛诊断和贝叶斯回归实战代码可直接复制运行。2. MCMC的马尔科夫链机制平稳分布、细致平衡与混合效率2.1 从“求不出积分”到“抽得出样本”蒙特卡洛在数据增强里的角色后验推断的通用形式是算后验期望 E[f(θ)|y] ∫f(θ)π(θ|y)dθ其中π(θ|y) ∝ L(y|θ)π(θ)。分子好写分母是那个归一化常数∫L(y|θ)π(θ)dθ维度一高就积不出来。蒙特卡洛的想法是从π(θ|y)直接抽取大量样本均值收敛到上述期望。问题在于“直接抽”在很多分布上做不到——逆变换要求可积的CDF拒绝采样在高维空间接受率指数级下降重要性采样在重尾时方差爆炸。MCMC换了个思路构造一条转移链让它在长期运行中“生成”目标分布的样本不需要求归一化常数也不需要直接采样。这条链的关键不在样本是否独立而在转移核P(x→x′)能使目标分布π在每一步之后都保持不变也就是不动点条件π πP。工程上更常用的是一个充分条件叫detailed balance细致平衡π(x)P(x→x′) π(x′)P(x′→x)满足它π一定是该链的平稳分布。MH算法的全部工作就是设计转移核让这个条件成立给定当前x先从建议分布q(x′|x)抽候选再以接受率决定是否移动。2.2 转移核与平稳分布MH如何满足detailed balanceMH的接受率形式是 α(x, x′) min(1, π(x′)q(x|x′) / (π(x)q(x′|x)))。如果提议分布对称比如高斯随机游走提议x′ x N(0, σ²)q(x|x′) q(x′|x)可以约掉接受率退化为min(1, π(x′)/π(x))这就是Metropolis采样。这个简化是matlab仿真里最常用的结构写目标分布时只需要定义联合密度而不用管归一化常数因为比值里常数消掉了。这正是MCMC最大的工程红利——很多应用问题里写得出似然×先验但归一化常数根本没闭式。接受率高低成了转移核设计里的核心变量。接受率过高意味着新样本离当前点很近链在局部徘徊混合慢过低意味着每步都在拒绝链原地踏步。高维情况下理论最优接受率在0.234附近Roberts等人对高斯目标的渐近结果一维目标大约0.44。仿真中不必精确逼近这个数只需监控接受率落在0.2~0.5区间偏离太多就调整提议方差。混合效率直接决定后面看到的自相关强度。2.2.1 接受率公式的对称性简化与数值风险细致平衡是局部条件它确保任意一步转移中从x到x′的期望流等于反向流。把两侧对x积分会发现π在转移作用下不变。这个证明在MATLAB中不必显式写出但它决定了接受率的结构分子分母必须包含前向概率π(x′)q(x′|x)和反向概率π(x)q(x|x′)二者比值失衡时通过随机接受修正。很多仿真发散的问题追根溯源是写接受率时把q项当成常量约掉破坏了详细平衡条件。2.3 burn-in、混合与自相关样本不是独立同分布的从任意初始值出发链需要一段“预热”才能进入平稳区域这一段称为burn-in样本要丢弃。进入平稳后相邻样本高度相关——这是马尔科夫链的固有属性不是bug。相关越强链的有效样本量越小。量化上使用自相关长度τ有效样本量ESS N/ττ越大信息越少。仿真时用trace plot横轴迭代次数纵轴样本值直接看混合状态理想情况下链在两个模式之间频繁穿梭而不是长时间困在一侧。三张图画完MCMC仿真才算真的“可见”了。3. 用MATLAB从零写Metropolis-Hastings采样器对准双峰混合高斯3.1 目标分布与提议分布的选型为什么用随机游走高斯提议双峰混合高斯是MCMC最常用的试验场分布有两个峰值如果链只在一个峰附近徘徊trace plot一眼就能看出来。设定目标分布为π(x) 0.6·N(x; -2, 0.8²) 0.4·N(x; 3, 1.2²)。这个分布没有解析采样方法但密度表达式简单。提议分布采用对称高斯随机游走x′ x σ·randn()。对称性的好处是接受率里没有q项实现简单缺点是高维时固定方差难以匹配不同维度尺度后面自适应MCMC章节会解决。在确定性仿真场景下σ的取值就是最重要的调参旋钮。先定义对数密度函数。这里有一个资深用户常踩的坑不要直接用exp相加再取log两个高斯分量强度差异大时exp(w)会下溢为0log变成-Inf。应该在log域计算用log-sum-exp技巧function logp log_target(x) % 双峰高斯混合的对数密度返回标量 % 每个分量的log权重 log高斯核归一化常数在MH比值中抵消 logw zeros(size(x, 1), 2); % 第一个分量均值-2标准差0.8权重0.6 logw(:, 1) log(0.6) - 0.5 * ((x 2) / 0.8).^2; % 第二个分量均值3标准差1.2权重0.4 logw(:, 2) log(0.4) - 0.5 * ((x - 3) / 1.2).^2; % log-sum-exp先减最大值再取exp防止指数下溢 m max(logw, [], 2); logp m log(sum(exp(logw - m), 2)); end这段函数的关键在于log-sum-exp把公共最大值m提出来exp(logw - m)的值域被控制在(0,1]log(sum)不会爆炸。均值、标准差、权重是目标分布的物理意义参数而- log(0.8·sqrt(2π))这类常数可以放心丢掉因为MH接受率是比值。3.2 主循环代码与接受率log域计算避免下溢MH主循环如下。以x00为初始状态前5000次作为burn-in丢弃后20000次保留% 主循环参数 N 20000; % 采样总数 burnin 5000; % 预热步数 sigma 1.5; % 提议标准差随机游走步长 x 0; % 初始状态 samples zeros(N, 1); acc 0; % 接受计数器 for t 1:(N burnin) % 高斯随机游走提议对称分布q项在MH中抵消 x_star x sigma * randn(); % log域接受率接受概率 exp(logp(x_star) - logp(x)) log_alpha log_target(x_star) - log_target(x); if log(rand()) log_alpha x x_star; acc acc 1; end % 丢弃burn-in阶段其余样本保存 if t burnin samples(t - burnin) x; end end acceptance_rate acc / (N burnin); fprintf(接受率: %.3f\n, acceptance_rate);主线是“提议-接受/拒绝-记录”。log(rand()) log_alpha把接受概率转化为一次均匀随机数与阈值的比较避免了对alpha再取exp的数值风险。计数器放在if内部分母包含burnin因为预热阶段的接受率同样反映提议方差是否合适。运行结束后接受率约在0.3附近说明步长合理低于0.1要减小sigma高于0.6则增大sigma。3.2.1 提议方差σ怎么设三个档位与接受率对照sigma的选择直接决定采样质量经验上分成三个档位sigma取值接受率区间自相关表现典型问题0.1~0.30.7相邻样本高度相关ESS小链移动太慢burn-in极长0.8~2.00.2~0.5适中trace快速穿梭最理想区间5以上0.05链长时间停在原地大量拒绝仿真发散感强提示自相关强不是采样失败而是有效样本量不足的信号。此时优先调sigma而不是盲目增加总迭代数。3.3 采样结果的三个验证图trace、直方图、自相关画验证图的代码集中在一次figure里subplot(2, 2, 1); plot(1:N, samples); title(Trace Plot); xlabel(迭代次数); ylabel(x); subplot(2, 2, 2); histogram(samples, 60, Normalization, pdf); hold on; % 把真实混合分布画上去做对比 xs linspace(-6, 7, 400); true_pdf 0.6 * normpdf(xs, -2, 0.8) 0.4 * normpdf(xs, 3, 1.2); plot(xs, true_pdf, r-, LineWidth, 2); legend(MCMC样本, 真实密度); subplot(2, 2, 3); autocorr(samples, 50); % 自相关图横轴lag纵轴自相关 fprintf(样本均值: %.3f\n, mean(samples)); fprintf(样本方差: %.3f\n, var(samples));3.3.1 三个图分别回答什么问题trace plot回答“链收敛了吗”。收敛的链纵向抖动幅度稳定两个峰之间周期性跨越若长时间只在一个峰附近要么是burn-in不够要么是sigma太小。直方图回答“样本分布对吗”。MCMC样本的直方图与真实密度红色曲线贴合说明抽样系统正确贴合但形状粗糙说明样本数不够或自相关太高。自相关图回答“信息量真的够吗”。自相关从lag0的1.0衰减到0的速度越快越好若滞后30仍然高于0.5有效样本量会显著低于名义样本量。操作视频的演示里通常把这三张图并列放在同一个figure窗口运行结束一眼确认收敛状态。下一步顺着变异方向走当变量维度增多靠手调sigma不现实就要引入Gibbs采样和自动协方差调节。4. 从Gibbs到自适应MCMC样本效率、R-hat与Geweke诊断4.1 Gibbs采样用满条件分布更新每一个分量当目标分布是多维联合分布而每个分量的满条件分布有闭型时Gibbs采样比MH更高效——不需要调节提议方差也没有拒绝步。做法是把θ拆成(θ1, …, θd)每一轮依次从p(θi|θ_{-i}, y)采样替换当前值。这个循环保持每个条件分布为平稳分布组合起来就采样了联合后验。工程上最常见的场景是共轭先验下的贝叶斯模型正态-正态、正态-逆伽马组合都有闭式条件分布。MATLAB中实现Gibbs的骨架如下% Gibbs采样骨架三变量目标条件分布假设均有闭型 theta zeros(3, 1); n_iter 10000; samples_gibbs zeros(n_iter, 3); % 伪代码示意实际条件分布需按模型代入 for iter 1:n_iter % 依次更新每个分量直接以更新的值作为后续分量条件 theta(1) conditional_sample_1(theta(2), theta(3)); theta(2) conditional_sample_2(theta(1), theta(3)); theta(3) conditional_sample_3(theta(1), theta(2)); samples_gibbs(iter, :) theta; end更新顺序有讲究第1个分量总是使用最新的θ2、θ3第2个分量使用刚更新的θ1和旧的θ3。这种“始终用新值”的扫描方式称为systematic scan是MATLAB中条件更新矩阵运算的标准写法能加快混合速度。4.2 自适应提议2.38²/d缩放与协方差学习期Gibbs要求条件分布可采样但很多后验并不满足。回到MH框架在多维情况下把单变量随机游走升级为多元高斯提议x′ x N(0, σ²·Σ)其中Σ是目标分布的协方差估计σ是缩放因子。常用的经验法则是σ 2.38 / sqrt(d)d为目标维度这个缩放因子在高斯目标下渐近最优。MATLAB里通过chol分解生成多元高斯样本% 自适应提议的协方差学习过程 d 2; Sigma_adapt eye(d) * 0.02; % 初始小协方差 L_adapt chol(Sigma_adapt, lower); N_total 8000; learn_period 3000; % 学习期之后提议固定 x zeros(d, 1); samples_adapt zeros(N_total, d); for t 1:N_total % 用当前协方差矩阵的Cholesky因子生成多元高斯增量 increment L_adapt * randn(d, 1); x_star x 2.38 / sqrt(d) * increment; log_alpha log_target_vec(x_star) - log_target_vec(x); if log(rand()) log_alpha x x_star; end samples_adapt(t, :) x; % 仅在学习期内更新协方差估计学习期结束后固定提议 if t learn_period mod(t, 50) 0 Sigma_adapt cov(samples_adapt(max(1, t-500):t, :)); Sigma_adapt Sigma_adapt 1e-8 * eye(d); % 正则化防奇异 L_adapt chol(Sigma_adapt, lower); end end这里两个必调参数学习期learn_period和正则项系数1e-8。学习期必须小于总迭代数结束后提议协方差固定否则链的非马尔科夫性破坏收敛理论。正则化防止协方差矩阵在样本少时奇异导致chol失败。注意cov()接收的是过去500个样本形成的窗口滑动窗口比全历史更能捕捉局部协方差结构。提示自适应MCMC只是在学习期用过去的样本来“调参”正式采样仍要求提议不变。仿真中若看到学习期之后接受率骤降多半是learn_period设置过短协方差还没收敛就锁定了。4.3 收敛诊断三件套R-hat、有效样本量与Geweke检验把MCMC当黑箱用时收敛诊断是保险。R-hat需要至少两条独立链比较链间方差与链内方差function Rhat compute_rhat(chains) % chains: M条链 × N次迭代矩阵 [M, N] size(chains); % 链均值的方差 → 链间方差B chain_means mean(chains, 2); overall_mean mean(chain_means); B N / (M - 1) * sum((chain_means - overall_mean).^2); % 每条链内方差的平均 → 链内方差W W mean(var(chains, 0, 2)); % 总方差估计混合公式 var_est (N - 1) / N * W B / N; Rhat sqrt(var_est / W); end计算时注意var(x, 0, 2)的维度参数0表示除以N-12表示按行计算。R-hat越接近1越好。有效样本量ESS针对单条链的自相关结构function ess effective_sample_size(x) % 通过自相关估计自相关长度tau再得到ESS n numel(x); xc xcorr(x - mean(x), coeff); rho xc(n:end); % 保留非负lag一侧 % 截断自相关尾部噪声对长lag求和影响大 max_lag min(1000, n - 1); tau 1 2 * sum(rho(2:max_lag 1)); ess n / tau; end自相关尾部噪声会让tau被高估max_lag上限正是为了压制这一点。Geweke检验取前10%与后50%的样本均值比较构造的z统计量落在[-1.96, 1.96]内表示未检测到收敛异常。三者的用法对照指标阈值/参考值结果解读R-hat1.1多链一致可认为已收敛ESS400后验均值估计稳定Geweke z|z| 1.96前后段均值无显著差异仿真发散时先看这三项再决定是加长学习期、增大迭代数还是回头检查接受率和提议方差。5. 实战闭环用MCMC拟合贝叶斯线性回归并输出后验可信区间5.1 模型设定与条件后验推导把前面的技术落到最常见的回归场景。模型为y Xβ εε ∼ N(0, σ²I)。配共轭先验β ∼ N(β0, Λ0⁻¹)σ² ∼ IG(a0, b0)。该设定下条件后验有闭型β|σ², y ∼ N( (XᵀX/σ² Λ0)⁻¹(Xᵀy/σ² Λ0β0), (XᵀX/σ² Λ0)⁻¹ )σ²|β, y ∼ IG( a0 n/2, b0 (y - Xβ)ᵀ(y - Xβ)/2 )。第二个条件分布只依赖残差这是Gibbs能闭环的关键。生成仿真数据rng(42); n 200; X [ones(n, 1), randn(n, 1), randn(n, 1)]; % 截距 两个自变量 beta_true [0.5; -1.2; 0.8]; y X * beta_true 0.5 * randn(n, 1);5.2 Gibbs循环的MATLAB实现按条件更新β和σ²% 共轭先验参数 d size(X, 2); beta0 zeros(d, 1); Lambda0 eye(d) * 0.01; a0 0.01; b0 0.01; % 初始值 beta beta0; sigma2 1.0; n_iter 12000; burnin 2000; beta_samples zeros(n_iter - burnin, d); sigma2_samples zeros(n_iter - burnin, 1); for iter 1:n_iter % 1. 更新beta多元正态协方差矩阵来自XX/σ² Λ的逆 Xbar X * X / sigma2 Lambda0; beta_mean Xbar \ (X * y / sigma2 Lambda0 * beta0); Sigma_beta inv(Xbar); Lb chol(Sigma_beta, lower); beta beta_mean Lb * randn(d, 1); % 2. 更新sigma²逆伽马注意gamrnd第二个参数是scale不是rate resid y - X * beta; a_post a0 n / 2; b_post b0 0.5 * (resid * resid); sigma2 1 / gamrnd(a_post, 1 / b_post); if iter burnin beta_samples(iter - burnin, :) beta; sigma2_samples(iter - burnin) sigma2; end end两个常见的坑都在这里gamrnd(a, b)第二个参数是scale逆伽马的rate参数需要取倒数所以写1/b_postbeta_mean用矩阵左除Xbar \ (...)而不是inv(Xbar)*...数值稳定性更好尤其当Xbar接近病态时。5.3 输出后验结果并与OLS对比采样结束后直接计算后验统计量beta_mean_mcmc mean(beta_samples, 1); beta_ci quantile(beta_samples, [0.025, 0.975], 1); beta_ols (X * X) \ (X * y); fprintf(OLS: %.3f %.3f %.3f\n, beta_ols); fprintf(MCMC: %.3f %.3f %.3f\n, beta_mean_mcmc); fprintf(95%% CI:\n); disp(beta_ci);输出对照关系大致如下表实际运行会有小数位波动真实系数应落在CI内参数OLS估计MCMC后验均值95% CIβ0约0.47约0.48[0.28, 0.69]β1约-1.22约-1.22[-1.29, -1.15]β2约0.83约0.82[0.74, 0.89]这里的核心收益是“输出完整后验而不只是点估计”OLS给出点和标准误MCMC给出分布可以随时改置信水平、计算P(β1 0)这类问题。收尾时跑一下effective_sample_size(beta_samples(:,2))把ESS卡在400以上再汇报结果若低于阈值就增加迭代数而不是缩短burn-in。整套流程对照操作视频逐帧确认burn-in段恰好对应trace plot从初始值摆动到平稳区域的时刻协方差学习期结束对应接受率由震荡转平缓的转折点这两个视觉标志对齐了你手上这套基于MCMC马尔科夫-蒙特卡洛抽样matlab仿真才算真正跑通。本文还有配套的精品资源点击获取