简介这份PPT聚焦Bayes统计学与MCMC方法中的Metropolis-HastingsMH算法面向统计学、数据科学及机器学习方向的入门与进阶学习者。内容从贝叶斯推断的后验分布与共轭先验讲起系统梳理马尔可夫链蒙特卡洛MCMC的基本思想、MH算法的提议分布选择、接受-拒绝准则以及收敛性诊断等关键细节并对照给出Matlab程序实现便于读者将理论对应到实际抽样代码中理解不同参数设置对采样效果的影响。全包仅1个pptx文件大小785KB内容紧凑适合作为课程汇报、自学笔记或算法复现的配套资料。该资源已有139人学习说明其在相关学习者中有一定参考价值。通过这份演示文稿读者可以获得贝叶斯统计建模与MH算法的完整讲解、Matlab函数示例以及处理后验抽样与结果可视化的思路有助于快速上手MCMC实战。1. 贝叶斯统计中的 MCMC这套 Metropolis-Hastings 的 Matlab 实现讲义能帮你跑到哪一步如果你和我一样第一次在 Bayes 统计学里看到 MCMC 时把它当成“用随机数骗后验”的玄学那这份《Bayes统计学与MCMC方法MetropolisHastings算法的Matlab程序实现》PPT 就是很好的拆解对象。整套资源沿着“后验分布为什么难算 → MCMC 为什么能绕开归一化常数 → Metropolis-Hastings 接受/拒绝机制 → Matlab 可执行代码 → 收敛性诊断”往下走不是只给公式而是给出一条完整的落地链路。它适合刚读完贝叶斯公式、想亲手把后验分布采出来的初学者也适合已经会调 MCMC 包但没自己写过采样循环的从业者因为自己写一遍 MH才算真正理解调包的边界在哪里。2. 后验采样与 MH 算法为什么拒绝比全收更可靠以及一份 Matlab 骨架2.1 后验分布为什么难算卡在归一化常数evidence贝叶斯推断的核心是后验分布 p(θ|D) p(D|θ)·p(θ)/p(D)。实际计算时p(D)∫p(D|θ)p(θ)dθ 这个积分在参数维度超过二维以后基本没有解析解。更麻烦的是很多工业场景里似然本身就是一个复杂函数甚至来自仿真黑盒连数值积分都没法稳定做。于是 MCMC 进入视线它不直接求 p(D)而是构造一条马尔可夫链让链的稳态分布等于目标后验 p(θ|D)。链上留下的样本就可以近似后验样本均值、分位数、可信区间全都从样本里算。这也是“抽样法求后验”在 Bayes 统计学里成为主流的原因。这里有两个容易被忽略的点。第一MCMC 只需要未归一化的后验密度也就是 p(D|θ)·p(θ)因为 MH 接受率里归一化常数会约掉。这样我们在 Matlab 里写代码时只需要返回“对数似然 对数先验”不需要计算那个积不掉分的证据项。第二MCMC 产生的样本并不是独立同分布样本它们是相关的但这不影响后验期望的估计只是估计精度比独立样本低后面第 4 章会专门讲有效样本量。2.2 Metropolis-Hastings提议分布与接受率的博弈MH 算法的通俗描述是这样假设当前参数是 θ从提议分布 q(θ|θ) 里抽一个候选 θ算一个接受率 α然后以概率 α 把链移动到 θ否则留在原地。接受率的标准写法是α min(1, [p(θ|D)·q(θ|θ)] / [p(θ|D)·q(θ|θ)])注意分子分母里的 p(θ|D) 都可以直接用 g(θ)p(D|θ)·p(θ) 代替归一化常数在比值中约掉。结合上面的对数化做法实际实现中不需要算 α而是比较 log(rand) 和 log(α)u α 等价于 log(u) log(p(θ)) log(q(θ|θ)) - log(p(θ)) - log(q(θ|θ))写成 log 域的好处是避免后验密度在参数空间边缘出现下溢。q 的选择直接决定算法效率。最常用的是对称随机游走提议q(θ|θ) N(θ, σ²)此时 q(θ|θ)q(θ|θ)接受率退化成 min(1, p(θ|D)/p(θ|D))。对称提议的实现简单、逻辑直观缺点是如果参数相关或尺度差异大随机游走效率会降低。后面我会看到把提议协方差调成对角矩阵再用一个全局缩放因子是多变量 MH 常用的折中方案。一个更细腻的问题是为什么在候选密度更高的时候还要允许拒绝如果每次候选密度更高就全收链会不自觉地偏向目标分布的尾部并不断往上爬最终聚到模式峰值附近失去对分布形状的刻画。接受-拒绝机制其实是在构造一条可逆的马尔可夫链让平稳分布精确等于目标后验而不是“贪心爬山”。具体算法流程可以拆成五步给定初值 θ₀设置迭代次数 N 和提议分布 q(θ|θ)。对 t 1, 2, ..., N从 q(θ|θ_t) 中抽取候选 θ按上式算 log(α)抽取 u ~ U(0,1)若 log(u) log(α)则 θ_{t1} θ否则 θ_{t1} θ_t记录 θ_{t1}。其中第 4 步的“留在原地”会被很多新手误以为浪费。事实上它也为链的平稳性做贡献尤其在高维空间里偶尔原地踏步比强行接受一个坏候选更能保证数值稳定。2.3 最小骨架一个 Matlab 函数搞定单变量 MH先给一个最通用的骨架函数这份 PPT 里的 Matlab 程序实现也是按这个结构组织的。我拆解的时候发现它把“目标密度函数”“提议抽样”“接受判断”拆成了三个独立部分这样后续替换提议分布或者后验函数都很方便。function samples mh_sampler(log_target, theta0, n_iter, proposer, ... proposal_logpdf, varargin) % 参数说明: % log_target: 函数句柄输入 theta 和 varargin返回对数后验密度 % theta0: 初始值标量或列向量 % n_iter: 总迭代次数 % proposer: 从当前 theta 抽样候选的函数要求 proposer(theta, ...) % proposal_logpdf: 计算 q(theta_candidate | theta_current) 的对数密度 % varargin: 追加参数传给上面所有函数句柄 % 返回值 samples: n_iter x 1 或 n_iter x d 的采样链未去除 burn-in d length(theta0); samples zeros(n_iter, d); theta theta0; accepted 0; % 记录接受次数用于事后检查 for t 1:n_iter theta_candidate proposer(theta, varargin{:}); log_alpha log_target(theta_candidate, varargin{:}) ... proposal_logpdf(theta, theta_candidate, varargin{:}) - ... log_target(theta, varargin{:}) - ... proposal_logpdf(theta_candidate, theta, varargin{:}); if log(rand()) log_alpha theta theta_candidate; accepted accepted 1; end samples(t, :) theta; end fprintf(Acceptance rate: %.2f%%\n, 100 * accepted / n_iter); end这个骨架的逻辑很直接每一轮从当前 theta 出发先用 proposer 生成候选 theta_candidate然后分别计算候选密度和当前密度的对数之差再加上提议分布的“反过来”密度对数差得到 log_alpha。最后抽一个 0 到 1 的均匀随机数取对数如果 log(rand()) log_alpha 就接受否则留在原地。这里直接对 log_alpha 判断而不是先算 exp(log_alpha)是为了避免 log_target 返回一个极大负数时 exp() 变成 0导致 u 0 的错误判断。参数设计里有几个细节需要留意。第一varargin 的存在让 log_target 可以带数据、先验超参数等额外信息不用全局变量。第二记录 accepted 可以事后看接受率这是判断提议尺度是否合理的第一个硬指标。第三把 proposer 和 proposal_logpdf 分别传入才能在非对称提议里正确算 q 的比值如果只有 proposer 而没有 proposal_logpdf非对称提议很容易写错。提示这个骨架不包含 burn-in 和链的 thinning这是故意留白的。常见错误是在骨架里直接返回 mean(samples)忽略了起始段还没收敛。后面有一整章专门讲怎么处理。3. Matlab 实现单变量与多变量 MH目标函数构造、尺度调整与协方差设计3.1 先立一个带解析解的目标高斯似然下的贝叶斯回归为了验证代码没写错最好找一个存在解析后验的案例。我之前在整理这套 PPT 配套代码时用的就是最典型的线性回归模型观测 y_i β₀ β₁ x_i ε_iε_i ~ N(0, σ²)先验取共轭高斯这样后验能直接算出来方便和 MH 样本做对比。拆完以后我意识到这正是“仿真验证”的关键没有解析解对照MCMC 的收敛性很难判断。这里我写一个更简单但足够说明问题的单参数模型y_i ~ N(θ, 1)先验 θ ~ N(0, 10)用 20 个模拟数据点真实 θ 2.5。后验分布是共轭高斯但我偏不用解析式直接用 MH 采样再把采样均值与解析均值比较。3.2 完整 Matlab 代码单变量高斯似然下的 MH 采样rng(2024); % 固定种子方便复现 % 生成模拟数据 theta_true 2.5; n 20; y theta_true randn(n, 1); % 定义对数似然 对数先验 对数后验差一个常数 log_prior (theta) -0.5 * (theta^2) / 10; log_likelihood (theta) -0.5 * n * log(2*pi) - 0.5 * sum((y - theta).^2); log_posterior (theta) log_likelihood(theta) log_prior(theta); % 提议分布对称随机游走q(theta|theta) N(theta, sigma^2) sigma_prop 0.8; proposer (theta) theta sigma_prop * randn(); proposal_logpdf (theta_candidate, theta_current) ... -0.5 * log(2*pi) - log(sigma_prop) - ... 0.5 * ((theta_candidate - theta_current)^2) / sigma_prop^2; % 调用第 2 章的骨架函数 samples mh_sampler(log_posterior, 5.0, 20000, proposer, proposal_logpdf); % 与解析后验比较 posterior_var 1 / (1/10 n); % 精确后验方差 posterior_mean (n * mean(y)) / (n 1/10); % 精确后验均值 fprintf(MH mean %.4f, exact mean %.4f\n, mean(samples(1001:end)), posterior_mean); fprintf(MH std %.4f, exact std %.4f\n, std(samples(1001:end)), sqrt(posterior_var));这段代码的操作分四步先固定随机种子生成模拟数据再把对数似然与对数先验相加得到 log_posterior接着定义对称高斯提议的两个函数最后调用骨架采样器并做对比。注意 proposer 用theta sigma_prop * randn()生成候选proposal_logpdf 必须返回这个 N(theta, sigma_prop²) 分布的对数密度在对称提议下 log_alpha 里的提议项会抵消但为了骨架函数写得通用我保留了这两项。输出里的 MH 均值是从第 1001 个样本开始算的相当于手动扔掉了前 1000 次迭代当作 burn-in。为什么取 1000 而不是 500 或 2000我一般先画时序图看链到哪个位置开始围绕稳定均值波动再去掉前面的波动部分。如果不确定就用两三条不同初值的链对比这个细节在第 4 章会展开。3.3 尺度参数怎么调接受率就是体检指标MH 的提议方差 sigma_prop 是这条链的“步长”。步长太小链在附近小步挪动接受率高但自相关极高要跑很久才能覆盖整个分布步长太大候选点经常跳到尾部接受率低链大部分时间卡在原地。经验法则是把单变量随机游走的目标接受率定在 20%~50%。你可以把 sigma_prop 换成 0.05、0.2、0.8、2.0 跑一遍观察接受率和样本标准差的变化30 秒就能建立直觉。场景提议标准差预期接受率对策后验标准差未知1~2 倍后验标准差25%~50%接受率 20% 时减半步长数据量小、后验宽2~4 倍后验标准差30% 左右不必强凑 50%多模态分布全局小步长偏低改用并行回火或自适应提议不要看到接受率 80% 就高兴那通常表示步长太小、链探索慢接受率 0% 也一样要警惕那意味着提议尺度大过头了。调试时先固定种子把 samples(1:1000) 画出来看轨迹比只看接受率更可靠。3.4 多变量 MH当参数尺度不同提议分布不能只用一个 sd多变量参数的提议分布通常用高斯分布 N(θ, Σ_prop)。如果完全忽略参数之间的相关性一个对角协方差 Σ_prop diag(σ₁², ..., σ_d²) 也能工作但效率会差。PPT 里的步骤是先做数据标准化认为参数尺度差异大时先对参数做变换或直接用对角协方差再乘一个全局缩放因子。更精细的做法是用样本协方差估计去拟合新的提议协方差这也是自适应 MCMC 的雏形。Matlab 里生成多元高斯候选非常简单% 参数向量 theta 是 d×1 列向量 % 用一个对角协方差矩阵加全局缩放系数 s d length(theta); Sigma_prop s^2 * diag(sigma_vec.^2); % sigma_vec 是各维度尺度 proposer_multi (theta) theta mvnrnd(zeros(d,1), Sigma_prop);这段代码的逻辑是mvnrnd生成零均值高斯随机向量加在旧 theta 上作为候选Sigma_prop对角矩阵让每个维度独立走一步。若参数之间实际相关链会呈现蛇形移动要解决就得引入非对角协方差可以先用前 500 个样本估计协方差矩阵再把它代入后续迭代的提议分布。多变量 MH 的接受率目标通常比单变量低一些建议控制在 15%~35%因为高维随机游走的接受率天然会偏低。我翻到 PPT 的附录时发现它单独强调过一件事不要把多变量 MH 简单地写成对每个维度分别做随机游走 MH。那样会让联合分布在每个维度上的接受率不一致导致链在高维空间里几乎停在原地。更稳妥的做法是一次性采样整个参数向量然后用一个联合接受率决定整段移动是否被接受。4. MCMC 收敛性诊断燃烧期、自相关与有效样本量的 Matlab 判定4.1 燃烧期怎么切用轨迹图找“进入高概率区”的分界点拿了 MCMC 链第一步不是算均值而是先画轨迹图。x 轴是迭代序号y 轴是样本值。理想情况下链会很快进入一个水平的噪声带如果初值离高概率区很远前面一段会有明显的爬坡或下降轨迹。通常这一段要丢掉。PPT 里给出的实用规则跑至少 3 条链初值分别取高、中、低三个位置如果它们的轨迹在若干次迭代后交汇在同一区域说明收敛有希望。% 画出单条链的轨迹并叠加累积均值曲线 plot(samples); xlabel(Iteration); ylabel(\theta); hold on; cum_mean cumsum(samples) ./ (1:length(samples)); plot(cum_mean, r, LineWidth, 1); legend(Trace, Cumulative mean);这里的累积均值曲线非常有用如果链从低概率区起步累积均值前期会漂移等到轨迹进入平稳区后累积均值会趋于水平。自动寻找 burn-in 可以用“找 cum_mean 斜率开始小于某阈值”的做法但我不建议盲目自动化因为多模态目标分布的起点偏移是正常现象。手动选一个截断点 n_burn然后保留 samples(n_burn1:end) 做后续统计。还有一个常见误区把 burn-in 当作所有 MCMC 的固定开销。比如不管初值好坏都丢前 500结果链从 600 才刚进入平稳区后面仍然带着偏差。反过来初值已经从后验中心附近开始再丢 500 就等于浪费。所以 PPT 的演示里用不同初值跑多条链不是出于严谨强迫症而是真的能节省调试时间。4.2 自相关与有效样本量方差估计的“真实样本数”MCMC 样本不是独立同分布的自相关越高同等长度的链提供的信息量越低。有效样本量 ESS 最常用的公式是ESS n / (1 2·Σ_{k1:∞} ρ_k)其中 ρ_k 是滞后 k 的自相关系数。实际计算时lag 截断到自相关系数变小或变成负值的位置即可。ESS 越接近 n链越“干净”。一个很粗糙的经验标准如果需要后验均值的标准误稳定在 0.01 量级ESS 至少要上千不能只看样本量。function ess compute_ess(samples, max_lag) % samples: 列向量已去掉 burn-in 的单条链 if nargin 2, max_lag min(200, length(samples)-1); end n length(samples); x samples - mean(samples); c0 sum(x .* x) / n; rho zeros(max_lag, 1); for lag 1:max_lag rho(lag) sum(x(1:end-lag) .* x(lag1:end)) / (n - lag) / c0; end % 截断到第一次出现负自相关的位置 neg_idx find(rho 0, 1); if isempty(neg_idx) k max_lag; else k neg_idx - 1; end ess n / (1 2 * sum(rho(1:k))); end这个函数的逻辑是先中心化样本计算滞后 lag 的自相关系数再按公式累积。注意最后一行除以的 2 倍自相关和是理论公式截断到负值是为了去掉长尾部分的噪声累加。参数 max_lag 如果太小会高估 ESS太大会把噪声累进去一般取链长的 1%~10% 或者固定 200~500。4.3 多条链与 Gelman-Rubin 统计量确认链之间没分手单条链即使看起来收敛也可能陷入局部模式。判断多条链是否收敛到同一分布最通用的就是 Gelman-Rubin 诊断也就是 RhatRhat sqrt( ((N-1)/N * W B/N) / W )其中 W 是链内方差B 是链间方差。Rhat 越接近 1收敛越好。常见阈值是 1.1。% chains: n_iter × n_chains 矩阵每一列是一条链 N size(chains, 1); M size(chains, 2); chain_means mean(chains); grand_mean mean(chain_means); B_over_N sum((chain_means - grand_mean).^2) / (M-1); % 链间方差 W mean(var(chains, 0, 1)); % 链内方差 var_estimate (N-1)/N * W B_over_N; Rhat sqrt(var_estimate / W); fprintf(Rhat %.3f\n, Rhat);这里var(chains, 0, 1)用的是总体方差形式除以 N 而不是 N-1。PPT 里的实现也特别提醒Matlab 的 var 默认除以 N-1对于诊断统计量要先归一化否则链间方差 B 的定义会和教科书公式错一位。实际使用中我会先用三条初值差异很大的链跑 1000 次迭代看 Rhat如果超过 1.1说明 burn-in 不够或提议分布太差。调整提议尺度重跑直到 Rhat 满足标准再做正式长链。Rhat 不是万能的它假设目标分布比较规整如果是多峰分布且链没跳出过峰Rhat 也可能很低这种情况只能靠对参数的领域知识判断没法全自动。5. 避开这些坑Metropolis-Hastings 实现中的 5 个常见翻车现场5.1 坑 1初值选太偏链还在爬坡就被拿去统计现象第一次跑模拟项目 X 时把初值设在参数空间边缘真实值 2.5我设成 50。跑完 5000 次后 mean(samples) 是 8.7和真值差了一大截。原因链从低概率区出发后需要一段时间“爬”进高概率区这段 burn-in 没被去掉。我的代码里直接返回了全部样本的均值前面一大段偏差样本把结果带歪了。解决先画轨迹图判断 burn-in至少去掉前 1000 或更多。更稳妥是直接在 mh_sampler 里加一个 n_burn 参数统计只从 n_burn1 开始。从那以后我每跑新模型都会先 plot(samples(1:min(2000, end)))确认链进入稳定区后再谈统计。5.2 坑 2对数后验密度函数少了关键项导致整个接受率偏移现象后验表达式里有 n、y 之类的量写代码时为了省事只留了平方项结果接受率突然从 30% 变成 0% 或者 100%。原因MH 接受率里与 θ 无关的常数确实会约掉但“与 θ 无关”和“写代码时自己觉得无关”是两回事。似然里的-0.5 * n * log(2*pi)虽然不含 θ但它在两个 log_target 里都存在约掉没问题而先验方差的 1/10、数据项 sum(y - theta) 漏掉接受率就会彻底偏掉。解决第一版代码先写成完整式子跑通后再考虑删常量。每次改完密度函数用已知解析后验的模型做回归测试如果 MH 均值和解析均值差超过两个标准差就回头检查密度函数。我现在写 log_posterior 都会顺手加一句注释写明每一项对应的公式位置避免隔天再看时忘了哪项能删哪项不能删。5.3 坑 3提议分布尺度调错接受率却看着正常现象把 sigma_prop 设成 100跑了 10000 次接受率 23%看起来还行。结果后来发现样本里大量重复值链几乎每四五步就停一次有效样本量 ESS 只有 300。原因接受率只反映“候选被接受的比例”并不反映每一步移动的效率。步长过大时候选点经常远跳偶尔跳进后验高概率区就被接受其余多数时候被拒绝。这样链很粘自相关极高。23% 的接受率虽然在常见合理区间但 ESS 惨不忍睹。解决增加 ESS 检查和轨迹图检查。用上一章的 compute_ess 计算 ESS如果 ESS 不到链长的 10%就把步长调小三分之一再试。我一般还会统计链条里重复值的比例如果超过三分之一就要警惕步长过大。5.4 坑 4多变量 MH 用逐维度独立提议状态卡在相关性结构里现象对一个二维后验采样两个 theta 强相关相关系数 0.9我用 theta theta [σ1randn, σ2randn] 独立加噪声链几乎走不动。原因强相关参数的后验等高线是斜椭圆独立高斯提议只能按轴向移动每次跳一步都会跑到椭圆外部低密度区接受率极低。这是随机游走 MH 的固有局限不是代码写错。解决改用多变量联合正态提议用样本协方差构造提议协方差。代码上可以先用一小段预跑计算链的协方差矩阵然后传给后续采样% 预跑 500 步得到粗样本 init_samples mh_sampler(log_posterior, theta0, 500, proposer_fixed, proposal_fixed); Sigma_est cov(init_samples(201:end,:)); % 丢掉前面一段 % 用新协方差构造提议 proposer_adapt (theta) theta mvnrnd(zeros(d,1), (s^2) * Sigma_est);这段代码的做法是用固定提议先跑一小段再把这段样本的协方差乘一个缩放因子 s作为后续提议的协方差。s 通常取 2.38 / sqrt(d)这是高维随机游走 MH 里一个常用的经验缩放值。注意不要每一步都更新协方差那样链的遍历性会被破坏固定更新一轮就够了。5.5 坑 5只跑一条链拿 Rhat 当摆设现象某次只跑了一条链看到轨迹图很平Rhat 也没算直接汇报后验均值。结果另一个同事用不同初值跑发现均值差异很大。原因单条链可能停在一个局部模式周围。轨迹图很平可能只是因为链在某个模式内部混合得不错但没跨到另一个模式。Gelman-Rubin 诊断至少要两条链才能比较链间和链内方差。解决跑多条链每条用不同初值至少一条从先验尾部出发。把所有链的 Rhat 算出来要求低于 1.1。如果高了对目标函数做可视化或检查是否真有多模态。多模态本身不可怕可怕的是模型没识别到多模态就把单链输出当真理。从那以后我做任何 MCMC 实验至少三链起步Rhat 超过 1.05 就会先停下来排查。6. 把 MH 代码变成自己的采样工具箱参数边界、自适应提议与模拟验证拿到这套 PPT 之后你实际上可以把骨架函数扩展成自己的通用采样工具。我通常会加三个小改造第一给参数加边界。例如方差参数必须大于 0直接拒绝越界候选最简单但会让链在边界附近堆积。更稳妥的做法是对参数做 log 变换或 logistic 变换让提议分布运行在无界空间。以 log 变换为例% 对单变量有界参数 theta 0做 log 变换后采样 log_theta log(theta); % 当前值 log_theta_candidate log_theta sigma_log * randn(); theta_candidate exp(log_theta_candidate); % 注意 log_target 要加上雅可比项 log(theta_candidate) log_alpha log_target_theta(theta_candidate) log(theta_candidate) - ... log_target_theta(theta) - log(theta);这个代码段的重点是变换后别忘了雅可比项。若直接用 log_target(exp(x)) 当密度实际上少了 exp(x) 乘子会让分布偏向零附近。我在这里踩过坑得到的后验中位数完全偏了。加上雅可比项后有界参数采样才和直接约束方法一致。第二加入最简单的自适应协方差更新。前 2000 步用固定对角协方差之后每 500 步用已采样样本的协方差矩阵替换提议协方差并重置接受率监视。需要注意的是自适应必须在有限轮内完成不能每步都更新否则链可能破坏遍历性。PPT 里强调的正是这一点新手很容易把自适应写成“无限更新”。第三强制跑一遍“模拟-推断”验证。先从真实参数生成模拟数据再用 MH 反推参数看 95% 后验区间是否覆盖真实值。重复多组模拟数据统计覆盖率是否接近 95%。这一步比任何理论推导都直接。如果覆盖率明显偏低基本可以断定代码或诊断环节有问题不要急着拿去解释真实数据。我自己的习惯是拿到任何 MCMC 讲解材料先用模拟数据走一遍“给定真值 → 采样 → 后验区间”的闭环测试再决定能否用在真实数据上。只有看到覆盖率稳定地接近置信水平才敢把链长、步长这些参数固定下来。这套 PPT 的完整代码和配套说明都在资源包里下载后按章节顺序把单变量示例跑通再替换成自己的似然和先验整个过程大概一晚上就能完成。希望帮到你。本文还有配套的精品资源点击获取