2026/9/13 11:38:02
FEATURED · 精选文章 
简介一份面向机器学习研究生与算法工程师的资源包紧密围绕变分贝叶斯、粒子滤波及边缘粒子滤波三大主题配套徐亦达老师的系统课件与可直接运行的MATLAB代码适合希望从理论推导过渡到实际建模的学习者。压缩包共含44个文件、体积64.48MB其中27份PDF课件系统讲解概率图模型与非参数方法12个M脚本展示变分贝叶斯、粒子滤波、吉布斯采样等经典算法另附PPT讲义、Jupyter Notebook笔记及说明文档。该资源已有648人学习使用。除了理论讲解还通过大量可运行脚本详细演示了重采样、KL散度最小化、动态系统状态估计等关键实现并包含行业讲座材料、说明文档与Notebook示例便于读者在动手调试中巩固算法理解并迁移到自己的研究或工程项目中。1. 把变分贝叶斯、粒子滤波、边缘粒子滤波放在同一张图里学才算真正入门概率推断学期末复习机器学习时很多人会在同一页课件上看到三个名字变分贝叶斯VB、粒子滤波PF、边缘粒子滤波RBPF。乍看都是“贝叶斯推断的近似手段”但它们的近似对象完全不同一个把后验分布限制在参数化族里做优化一个用一堆带权重的样本去逼近后验另一个是前两者的混合体把能算的部分解析掉、只对剩下的状态采样。理解这三者的差异比背下公式重要得多。做状态空间模型的人最容易体会到这件事。给定观测序列想求状态的后验分布但非线性的状态转移或观测方程会立刻让后验失去闭式解。变分法适合“后验长什么样大体知道只是算不动”粒子滤波适合“后验形状完全未知但能从模型里采样”边缘粒子滤波则踩在中间适用于问题里恰好存在一块线性高斯结构。这篇内容把这三条路分别推到最小可运行的 MATLAB 实现上再讲清楚参数怎么设、边界在哪。适合正在做贝叶斯相关课题、准备机器学习期末复习、或者要在 MATLAB 里跑通第一个滤波示例的人。2. 从后验近似看变分贝叶斯ELBO、均值场分解与MATLAB实现2.1 变分贝叶斯在近似什么用ELBO把推断变成优化变分贝叶斯不直接计算后验 p(z|x)而是找一个分布 q(z) 去逼近它最常用的度量是 KL 散度。直接最小化 KL(q‖p) 还是绕不开后验的归一化常数于是把目标等价改写成最大化证据下界ELBOELBO(q) E_q[log p(x,z)] - E_q[log q(z)]这个式子同时包含“让 q 尽量解释数据”和“让 q 别太复杂”两个压力。如果把 log p(x) 看成模型证据ELBO 就是它减去一个非负的 KL 项所以 log p(x) ≥ ELBO名字里的“下界”由此而来。变分推断做的事就是去最大化这个下界。实际操作里最常用的假设是均值场分解把 q(z) 拆成各个因子 q(z_j) 的乘积忽略变量间在近似分布中的相关性。这个假设当然有代价但换来的是坐标上升的闭式更新每次只更新一个因子固定其余因子。跟过李宏毅机器学习公开课或者翻过任意一本概率图模型教材的人都会对这个“固定其他、轮流优化”的流程眼熟。它和机器学习里的梯度下降是同一类迭代思路只不过更新方向不是梯度而是每个因子的最优分布形式。近似误差从哪来一个来自均值场假设本身一个来自变分族表达能力的限制。这正好对应机器学习数学理论里泛化误差界的讨论框架偏差来自假设族不够宽方差则和样本量相关。2.2 一个能跑的变分贝叶斯例子高斯均值与精度估计下面用一个最简单的模型看清均值场 VB 的迭代结构数据 y_i 独立同分布于正态分布均值 μ 和精度 τ 都未知对 μ 给高斯先验、对 τ 给 Gamma 先验。这不是玩具高斯混合模型做贝叶斯推断时每个分量的参数更新就是这样一小块。rng(0); y randn(200, 1) * 0.5 1.2; mu0 0; lam0 1e-3; a0 1; b0 1; muN mean(y); lamN 1 / var(y); aN a0 numel(y) / 2; bN b0; for it 1:200 % 固定 q(mu)更新 q(tau)形状参数已定只更新速率参数 E_tau aN / bN; % 固定 q(tau)更新 q(mu)高斯后验的精度和均值 lamN lam0 numel(y) * E_tau; muN (lam0 * mu0 E_tau * sum(y)) / lamN; % 重新计算 q(tau) 的速率参数用到 q(mu) 的二阶矩 E_mu_mu0_sq (muN - mu0)^2 1 / lamN; E_y_mu_sq sum((y - muN).^2) numel(y) / lamN; bN b0 0.5 * (lam0 * E_mu_mu0_sq E_y_mu_sq); end这段代码的每一行都要能对上推导。E_tau 是当前 Gamma 分布下 τ 的期望更新 q(μ) 时用 E_tau 替代 τ 进入高斯似然项得到新的精度 lamN 和均值 muN。更新 bN 时则反过来需要 q(μ) 的二阶矩 E[μ²] (muN-mu0)² 1/lamN。坐标上升的对称性在这里看得最清楚两边互相依赖交替推进。参数里最容易调错的是先验的尺度。lam0 表示对 μ 先验的信心取 1e-3 代表“数据说了算”取 1 或更大则会把 μ 拉向 mu0。a0 和 b0 控制 τ 的先验期望Gamma 分布的均值为 a0/b0想让先验弱就同时取小值。初始值给 muNmean(y)、lamN1/var(y) 会让迭代快速稳定但注意它来自数据的矩估计并非先验信息如果刻意想观察收敛路径可以改成任意数再跑一次对比。2.3 迭代停止与初值这两个设置决定能不能收敛变分迭代没有“学习率”这个旋钮但有两个等价物迭代上限和参数变化容差。常见做法是设最大 100~300 次同时记录每轮 muN 的绝对变化量当变化小于 1e-6 就提前退出。如果把迭代过程想成坐标上升在 ELBO 曲面上爬坡那么和机器学习中的梯度下降有同一个直觉初值离局部最优太远前几十轮会走一大段“冤枉路”但不影响最终结果真正的风险是均值场假设把多峰后验逼成了单峰近似这时迭代再久也不会回到真实后验的形状。判定收敛的正确指标不是 muN 本身而是 ELBO。每轮多算一行下界值既能确认迭代在上升也能在后面 6.2 节作为实现正确与否的判据。如果 ELBO 出现锯齿状下降通常不是迭代次序的问题而是某个因子更新里的期望算错了。另一个常见坑是忘记更新 bN 里的numel(y)/lamN这一项少它一项τ 的估计会系统性偏大因为模型少算了 q(μ) 本身的不确定性。这个细节正是“课件推导能跳过代码跳不过”的地方。3. 粒子滤波用重采样对抗权值退化盯住三个环节3.1 从重要性采样到有效样本数Neff粒子滤波的思路是后验算不出来但模型可以向前模拟。假设从建议分布 q(x_t | x_{t-1}, y_t) 里抽出一批粒子每个粒子配一个权重更新规则是w_t ∝ w_{t-1} · p(y_t|x_t) · p(x_t|x_{t-1}) / q(x_t|x_{t-1}, y_t)。如果把状态转移先验当作建议分布更新就退化成“预测一步、按似然加权”。问题在于每轮只乘不除权重会迅速集中在极少数粒子上这一现象叫权值退化。衡量退化程度的量是有效样本数Neff 1 / sum(w.^2)取值为 1 到 N 之间。Neff 接近 N说明权重均匀Neff 掉到 N/2 以下说明大部分粒子已经不重要了继续推也只是浪费算力。重采样做的就是从当前粒子集合里按权重重新抽 N 个把权重重置均匀让资源重新分配到高概率区域。所以一个粒子滤波实现本质就是“预测、加权、重采样”三个环节的循环。3.2 SIR粒子滤波的MATLAB最小实现这里给出一个完整可跑的 SIR 粒子滤波函数状态和观测都是一维标量方便对着公式逐行查。function [x_est, samples] sir_pf(y, f, g, Q, R, N) % y: 观测序列行向量 % f: 状态转移函数句柄如 (x) 0.9*x % g: 观测函数句柄如 (x) x.^2/20 % Q: 过程噪声方差R: 观测噪声方差N: 粒子数 T length(y); x randn(1, N) * sqrt(Q); w ones(1, N) / N; samples zeros(T, N); x_est zeros(1, T); for t 1:T % 预测从状态转移先验中采样 x f(x) sqrt(Q) * randn(1, N); % 加权用观测似然更新权重 w w .* normpdf(y(t) - g(x), 0, sqrt(R)); w w / sum(w); % 重采样判断 Neff 1 / sum(w.^2); if Neff 0.5 * N idx systematic_resample(w); x x(idx); w ones(1, N) / N; end samples(t, :) x; x_est(t) sum(w .* x); end end function idx systematic_resample(w) % 系统重采样一次随机偏移产生 N 个等距采样点 N numel(w); edges ((0:N-1) rand) / N; edges(end) min(edges(end), 1 - eps); cw cumsum(w); idx zeros(N, 1); j 1; for i 1:N while cw(j) edges(i) j j 1; end idx(i) j; end end使用时先把观测数据准备好再定义转移和观测函数。normpdf来自统计工具箱如果电脑里没装手写exp(-0.5*((y-g(x))/sqrt(R)).^2) / sqrt(2*pi*R)效果完全一样。重采样阈值 0.5*N 是最常见的经验值取 0.3 会更省重采样次数取 0.7 会更激进样本多样性更好但计算量上升。注意每次重采样都会引入额外随机性所以必须保存随机种子这一点在 6.1 节会展开。系统重采样比多项重采样方差更小因为它只抽一个随机数其余 N-1 个点等距排在单位区间上。实现里的edges ((0:N-1) rand) / N就是那个“一次偏移加等距点”用累积分布函数的逆变换把每个区间映射成某个粒子的索引。最后的min(edges(end), 1-eps)是防边界越界最后一个点理论上不可能取到 1但浮点误差下还是兜底一下更稳。3.3 重采样策略与观测噪声很小的退化陷阱重采样选哪种影响的是粒子多样性。多项重采样每次独立抽 N 个序号方差大粒子重复率高系统重采样把 [0,1] 均匀切成 N 段每段取一个点重复率低。实际工程里系统重采样几乎总是更好的默认选项。另一个容易被忽略的是“重采样后要不要加微小抖动”如果需要维持粒子多样性可以给重复粒子叠一层极小幅度的过程噪声但这会引入额外偏差慎用。更隐蔽的坑来自观测噪声 R 的设置。R 设得比真实噪声小一个量级时似然函数变得非常尖权重在第一步就会全部集中到一两个粒子上Neff 瞬间掉到 2 以下。这时候重采样救不回来因为重采样只能复制“还不错的粒子”不能凭空生成“更好的粒子”。现象就是滤波轨迹出现成段的平直粒子全部相同后续预测只在同一条轨迹上加噪声。排查方法很简单把每轮的 Neff 画出来如果在前几步就枯竭先检查 R再检查观测函数是否写对了。状态转移噪声 Q 则相反Q 太小时粒子长期挤在一起多样性不足Q 太大时粒子发散靠似然拉回来需要更多样本。4. 边缘粒子滤波把可解析积分的状态边缘化减少粒子维度4.1 Rao-Blackwellization 为什么能压低方差边缘粒子滤波在英文里常叫 Rao-Blackwellized Particle Filter中文里的“边缘”指的是对部分状态做解析边缘化。做法是把状态拆成两块x_r 保留给粒子x_l 在给定 x_r 和观测的条件下可以解析求解通常是一个卡尔曼滤波问题。粒子只负责采样 x_rx_l 的条件后验用解析高斯分布表达。方差为什么变小可以看条件方差公式Var(X) E[Var(X|Y)] Var(E[X|Y])。用条件期望 E[X|Y] 去替代 X损失的只是条件方差那一项而这部分正是被解析滤波吸收掉了。对比纯粒子滤波RBPF 把一个高维采样问题拆成“低维采样 解析滤波”粒子维度一降Neff 的退化速度大幅放缓。这也是 RBPF 在目标跟踪、同步定位与地图构建里受欢迎的原因位置姿态用粒子地图特征或线性速度部分用卡尔曼各干各的。4.2 粒子采样与卡尔曼更新交替的RBPF主循环假设模型里 x_l 是标量且动态是线性的x_r 通过非线性函数影响 x_l 的演变。为了让结构清楚这里展示主循环的核心代码强调每个粒子都背着一个小卡尔曼滤波器。% RBPF 主循环骨架Np 个粒子P 保存每个粒子对 x_l 的方差 xl_pred zeros(Np, 1); P_pred zeros(Np, 1); for t 1:T % 第 1 步对 x_r 做粒子传播 xr fr(xr) sqrt(Qr) * randn(Np, 1); % 第 2 步对每个粒子做 x_l 的卡尔曼预测 for p 1:Np A a_func(xr(p)); % 线性系数随 x_r 变化 xl_pred(p) A * xl(p); P_pred(p) A^2 * P(p) Ql; end % 第 3 步观测更新与权重更新同时做 for p 1:Np innov y(t) - xl_pred(p); S P_pred(p) R; K P_pred(p) / S; xl(p) xl_pred(p) K * innov; P(p) (1 - K) * P_pred(p); % 权重用边缘似然即观测在高斯预测分布下的密度 w(p) w(p) * normpdf(innov, 0, sqrt(S)); end w w / sum(w); if 1 / sum(w.^2) 0.5 * Np idx systematic_resample(w); xr xr(idx); xl xl(idx); P P(idx); w ones(Np, 1) / Np; end end这段代码的关键在第三层循环中的权重更新权重乘的不是“给定 x_l 后验的似然”而是边际似然即只把 x_r 和观测之间的信息算进来。卡尔曼增益 K 负责修正 x_l而 x_r 的好坏完全由 innovation 的密度体现数值上正好是normpdf(innov, 0, sqrt(S))。如果这一步写成了用完整似然更新权重等于把 x_l 的信息重复计算了一遍滤波会表现得过于自信协方差估计偏小。重采样时最容易漏掉的是 P。xl 被重采样了但 P 还留在旧粒子上下一步的 P_pred 就会和 xl 不匹配卡尔曼增益失真。所以 P 必须跟着 xl 一起按 idx 重排。另一个结构上的经验是能解析的部分尽量拉大。x_r 里如果混进一个弱可观测的维度粒子数需求立刻上升RBPF 的优势就被削弱了。4.3 写RBPF最容易踩的三个代数坑第一状态增广的顺序。RBPF 里 x_r 在粒子中流动但 x_l 的协方差 P 是每个粒子各自维护的矩阵维度必须始终一致。很多人第一步写对了改动模型后 P 的维度忘了跟着改报错位置在卡尔曼增益计算处看起来像“除零”实际上是维度不匹配。建议在进入主循环前用assert(isequal(size(P), [Np, size(xl, 2)]))这类检查把所有数组维度卡死。第二权重更新用边缘似然而不是条件似然。数学上p(y_t | x_{r,1:t}, y_{1:t-1})是高斯分布均值为 xl_pred方差为 S也就是代码里的normpdf(innov, 0, sqrt(S))。如果这里误写成normpdf(y(t) - g([xr(p); xl(p)]), 0, sqrt(R))等于把 x_l 的解析结果也当成了采样变量MCMC 意义上没有错但方差会变大Neff 下降速度明显变快。判断方法是对比两种写法下的 Neff 曲线RBPF 的正确实现曲线应当更平滑。第三可观测性。x_l 的卡尔曼更新依赖观测方程里能看到它。如果某段时间观测对 x_l 完全无信息P 会不降反升权重也退化成只由 x_r 的预测密度决定。这种情况不是代码错而是模型本身在结构上不可观。取舍方法是把不可观的那部分状态挪到 x_r 里或者加一个惩罚先验。5. 调参顺序与选型粒子数、迭代次数与计算成本的取舍5.1 三种方法的适用边界先选结构再选参数对比项变分贝叶斯 (VB)粒子滤波 (PF)边缘粒子滤波 (RBPF)近似对象参数化分布 q(z)后验的经验分布x_r 经验分布 x_l 解析高斯对非线性支持弱依赖变分族灵活性强直接采样只对 x_r 支持非线性典型计算瓶颈迭代轮次 × 状态维度粒子数 × 步长粒子数 × (卡尔曼 积分)初始化敏感度高中中MATLAB 官方函数无通用 VB 函数部分版本工具集有 particleFilter需自行拼装选型逻辑从数据结构出发。如果模型里存在明显可解析的线性高斯块优先考虑 RBPF它的粒子维数最低同样的粒子数下精度最高。如果模型完全非线性且没有高斯结构PF 是兜底方案。VB 适合你要的不是逐时刻的状态轨迹而是参数的完整后验比如高斯混合模型的聚类中心后验它给的是一个分布形状而非样本序列。对做机器学习期末复习或要交课程作业的人来说最简单的判别方式是看关键变量是否连续可导连续且先验共轭用 VB强非线性且转移函数复杂用 PF两者混合用 RBPF。5.2 变分贝叶斯的迭代次数、容差与初值直觉VB 的迭代参数比粒子滤波少但更微妙。先设置一个最大迭代次数比如 300再设置 ELBO 变化量的容差 1e-4 或 1e-6。迭代次数不是越大越好因为均值场分解的偏差不会随迭代消失多跑几千轮只是逼近同一个目标。初值的影响集中在“对称性破缺”的模型上比如高斯混合模型如果初始把两个分量放在同一位置坐标上升可能永远不分开但这在简单的高斯均值模型上没有影响。还有一个从优化工具箱里带过来的习惯用目标函数的变化量做停止判据别直接用参数变化。ELBO 是目标muN 是参数。参数变化很小而 ELBO 还在爬坡说明进入了平坦区域ELBO 先升后降说明某一步期望算错。MATLAB 的优化工具箱在调试时能帮你验证凸问题下的收敛行为但 VB 的更新不是梯度下降不需要也不应该设置学习率。如果看到迭代震荡去检查 q(μ) 更新里是否用了最新的 bN而不是上一轮的旧值。5.3 粒子数与重采样阈值的工程经验值粒子数 N 的经典经验值一维状态 200~1000二维 1000~5000三维以上 5000~20000。这只是一个起点真正要看的指标是重采样频率和 Neff 最低值。如果 Neff 长期低于 0.2N先别急着加粒子而是调整过程噪声 Q 和观测噪声 R。增加 Q 会提高粒子多样性降低 R 会让粒子更集中多数情况下调噪声比重采样阈值有效。在学校实验室搭的机器学习服务器上跑大规模任务时粒子滤波的循环天然可并行每个粒子的预测和加权互不依赖。MATLAB 里可以把第一层 for 换成parfor但要注意重采样步骤必须在并行循环外统一做否则粒子索引会乱。RBPF 的每个粒子有独立的卡尔曼状态同样适合 parfor只要 P、xl、xr 都按切片传进去。阈值 0.5N 和 0.7N 之间的差别远小于噪声参数带来的差别不必在这上面过度纠结先把一套固定的设置跑通再观察 Neff 曲线决定往哪个方向调。6. 在MATLAB里验证这三种实现是否写对的三个技巧6.1 固定随机种子并保存中间状态让每次运行完全可复现粒子滤波和变分贝叶斯都吃随机数。粒子滤波的随机来自初始粒子和重采样VB 如果初始化用了randn同样不可复现。在脚本最上方加一行rng(2024)是最基本的但还不够调试时需要定位到“第几步的重采样出了问题”所以要保存每一步的随机数生成器状态。rng(0); seed_track cell(T, 1); for t 1:T seed_track{t} rng; % 粒子滤波或RBPF的预测与加权步骤 endrng不带参数返回值时取回当前生成器状态记录在 cell 里。下次出错后把生成器恢复成rng(seed_track{t})再单步运行粒子轨迹会完全复现。这个技巧对粒子滤波尤其重要因为重采样一旦混入不同机器或不同 MATLAB 版本产生的随机序列都不同定位到具体哪一步才能把问题从随机性中剥离开。6.2 用ELBO或对数似然曲线定位实现错误VB 实现最常见的错误是更新公式少了一项“‘次要项’曲线不会立即变红而是收敛到错误的 ELBO 值。正确写法是对每轮更新后的 q 分布重新计算完整 ELBO单独存一条曲线。如果 ELBO 最后稳定在某个值但没有全程单调递增说明至少一个因子更新前后没有提高下界问题多半出在该因子更新时用了过期参数。粒子滤波这边可以计算对数似然log p(y_{1:T})的估计公式为每一轮权值和的均值之积取对数。写一个监控变量log_lik log_lik log(mean(unnorm_w))如果它出现异常跳变通常是观测函数写反了或者观测噪声 R 数量级错了。这个方法比看状态轨迹更灵敏因为轨迹误差会被局部平滑掩盖而似然直接把每一步的预测质量摊开。6.3 让工具箱实现当对照互相对比状态估计轨迹MATLAB 较新版本里控制或系统辨识相关工具箱提供particleFilter对象输入which particleFilter可以确认本机是否可用。用它构建同样的模型跑一遍再与自己的 SIR 实现做轨迹对比。两者不应逐点完全一致因为随机性不同但状态估计曲线的趋势和波动幅度应当接近。如果自己的实现输出明显偏离最常见原因是重采样索引方向写错导致粒子被整体倒序重排。另一种对照是用unscentedKalmanFilter做参考线。无迹卡尔曼滤波对弱非线性有一阶精度粒子滤波在粒子数充足时更精确。两者之间的差距应该在合理范围内差距太小说明粒子数可能偏多或问题太线性差距太大优先查重采样和似然函数。最后一招是把过程噪声和观测噪声都设成极小值此时模型近似确定性粒子滤波的输出应该趋近于状态转移函数本身的轨迹。如果这个条件下还发散代码一定存在结构性错误与参数无关。本文还有配套的精品资源点击获取
变分贝叶斯、粒子滤波与边缘粒子滤波：从原理到MATLAB实现
简介一份面向机器学习研究生与算法工程师的资源包紧密围绕变分贝叶斯、粒子滤波及边缘粒子滤波三大主题配套徐亦达老师的系统课件与可直接运行的MATLAB代码适合希望从理论推导过渡到实际建模的学习者。压缩包共含44个文件、体积64.48MB其中27份PDF课件系统讲解概率图模型与非参数方法12个M脚本展示变分贝叶斯、粒子滤波、吉布斯采样等经典算法另附PPT讲义、Jupyter Notebook笔记及说明文档。该资源已有648人学习使用。除了理论讲解还通过大量可运行脚本详细演示了重采样、KL散度最小化、动态系统状态估计等关键实现并包含行业讲座材料、说明文档与Notebook示例便于读者在动手调试中巩固算法理解并迁移到自己的研究或工程项目中。1. 把变分贝叶斯、粒子滤波、边缘粒子滤波放在同一张图里学才算真正入门概率推断学期末复习机器学习时很多人会在同一页课件上看到三个名字变分贝叶斯VB、粒子滤波PF、边缘粒子滤波RBPF。乍看都是“贝叶斯推断的近似手段”但它们的近似对象完全不同一个把后验分布限制在参数化族里做优化一个用一堆带权重的样本去逼近后验另一个是前两者的混合体把能算的部分解析掉、只对剩下的状态采样。理解这三者的差异比背下公式重要得多。做状态空间模型的人最容易体会到这件事。给定观测序列想求状态的后验分布但非线性的状态转移或观测方程会立刻让后验失去闭式解。变分法适合“后验长什么样大体知道只是算不动”粒子滤波适合“后验形状完全未知但能从模型里采样”边缘粒子滤波则踩在中间适用于问题里恰好存在一块线性高斯结构。这篇内容把这三条路分别推到最小可运行的 MATLAB 实现上再讲清楚参数怎么设、边界在哪。适合正在做贝叶斯相关课题、准备机器学习期末复习、或者要在 MATLAB 里跑通第一个滤波示例的人。2. 从后验近似看变分贝叶斯ELBO、均值场分解与MATLAB实现2.1 变分贝叶斯在近似什么用ELBO把推断变成优化变分贝叶斯不直接计算后验 p(z|x)而是找一个分布 q(z) 去逼近它最常用的度量是 KL 散度。直接最小化 KL(q‖p) 还是绕不开后验的归一化常数于是把目标等价改写成最大化证据下界ELBOELBO(q) E_q[log p(x,z)] - E_q[log q(z)]这个式子同时包含“让 q 尽量解释数据”和“让 q 别太复杂”两个压力。如果把 log p(x) 看成模型证据ELBO 就是它减去一个非负的 KL 项所以 log p(x) ≥ ELBO名字里的“下界”由此而来。变分推断做的事就是去最大化这个下界。实际操作里最常用的假设是均值场分解把 q(z) 拆成各个因子 q(z_j) 的乘积忽略变量间在近似分布中的相关性。这个假设当然有代价但换来的是坐标上升的闭式更新每次只更新一个因子固定其余因子。跟过李宏毅机器学习公开课或者翻过任意一本概率图模型教材的人都会对这个“固定其他、轮流优化”的流程眼熟。它和机器学习里的梯度下降是同一类迭代思路只不过更新方向不是梯度而是每个因子的最优分布形式。近似误差从哪来一个来自均值场假设本身一个来自变分族表达能力的限制。这正好对应机器学习数学理论里泛化误差界的讨论框架偏差来自假设族不够宽方差则和样本量相关。2.2 一个能跑的变分贝叶斯例子高斯均值与精度估计下面用一个最简单的模型看清均值场 VB 的迭代结构数据 y_i 独立同分布于正态分布均值 μ 和精度 τ 都未知对 μ 给高斯先验、对 τ 给 Gamma 先验。这不是玩具高斯混合模型做贝叶斯推断时每个分量的参数更新就是这样一小块。rng(0); y randn(200, 1) * 0.5 1.2; mu0 0; lam0 1e-3; a0 1; b0 1; muN mean(y); lamN 1 / var(y); aN a0 numel(y) / 2; bN b0; for it 1:200 % 固定 q(mu)更新 q(tau)形状参数已定只更新速率参数 E_tau aN / bN; % 固定 q(tau)更新 q(mu)高斯后验的精度和均值 lamN lam0 numel(y) * E_tau; muN (lam0 * mu0 E_tau * sum(y)) / lamN; % 重新计算 q(tau) 的速率参数用到 q(mu) 的二阶矩 E_mu_mu0_sq (muN - mu0)^2 1 / lamN; E_y_mu_sq sum((y - muN).^2) numel(y) / lamN; bN b0 0.5 * (lam0 * E_mu_mu0_sq E_y_mu_sq); end这段代码的每一行都要能对上推导。E_tau 是当前 Gamma 分布下 τ 的期望更新 q(μ) 时用 E_tau 替代 τ 进入高斯似然项得到新的精度 lamN 和均值 muN。更新 bN 时则反过来需要 q(μ) 的二阶矩 E[μ²] (muN-mu0)² 1/lamN。坐标上升的对称性在这里看得最清楚两边互相依赖交替推进。参数里最容易调错的是先验的尺度。lam0 表示对 μ 先验的信心取 1e-3 代表“数据说了算”取 1 或更大则会把 μ 拉向 mu0。a0 和 b0 控制 τ 的先验期望Gamma 分布的均值为 a0/b0想让先验弱就同时取小值。初始值给 muNmean(y)、lamN1/var(y) 会让迭代快速稳定但注意它来自数据的矩估计并非先验信息如果刻意想观察收敛路径可以改成任意数再跑一次对比。2.3 迭代停止与初值这两个设置决定能不能收敛变分迭代没有“学习率”这个旋钮但有两个等价物迭代上限和参数变化容差。常见做法是设最大 100~300 次同时记录每轮 muN 的绝对变化量当变化小于 1e-6 就提前退出。如果把迭代过程想成坐标上升在 ELBO 曲面上爬坡那么和机器学习中的梯度下降有同一个直觉初值离局部最优太远前几十轮会走一大段“冤枉路”但不影响最终结果真正的风险是均值场假设把多峰后验逼成了单峰近似这时迭代再久也不会回到真实后验的形状。判定收敛的正确指标不是 muN 本身而是 ELBO。每轮多算一行下界值既能确认迭代在上升也能在后面 6.2 节作为实现正确与否的判据。如果 ELBO 出现锯齿状下降通常不是迭代次序的问题而是某个因子更新里的期望算错了。另一个常见坑是忘记更新 bN 里的numel(y)/lamN这一项少它一项τ 的估计会系统性偏大因为模型少算了 q(μ) 本身的不确定性。这个细节正是“课件推导能跳过代码跳不过”的地方。3. 粒子滤波用重采样对抗权值退化盯住三个环节3.1 从重要性采样到有效样本数Neff粒子滤波的思路是后验算不出来但模型可以向前模拟。假设从建议分布 q(x_t | x_{t-1}, y_t) 里抽出一批粒子每个粒子配一个权重更新规则是w_t ∝ w_{t-1} · p(y_t|x_t) · p(x_t|x_{t-1}) / q(x_t|x_{t-1}, y_t)。如果把状态转移先验当作建议分布更新就退化成“预测一步、按似然加权”。问题在于每轮只乘不除权重会迅速集中在极少数粒子上这一现象叫权值退化。衡量退化程度的量是有效样本数Neff 1 / sum(w.^2)取值为 1 到 N 之间。Neff 接近 N说明权重均匀Neff 掉到 N/2 以下说明大部分粒子已经不重要了继续推也只是浪费算力。重采样做的就是从当前粒子集合里按权重重新抽 N 个把权重重置均匀让资源重新分配到高概率区域。所以一个粒子滤波实现本质就是“预测、加权、重采样”三个环节的循环。3.2 SIR粒子滤波的MATLAB最小实现这里给出一个完整可跑的 SIR 粒子滤波函数状态和观测都是一维标量方便对着公式逐行查。function [x_est, samples] sir_pf(y, f, g, Q, R, N) % y: 观测序列行向量 % f: 状态转移函数句柄如 (x) 0.9*x % g: 观测函数句柄如 (x) x.^2/20 % Q: 过程噪声方差R: 观测噪声方差N: 粒子数 T length(y); x randn(1, N) * sqrt(Q); w ones(1, N) / N; samples zeros(T, N); x_est zeros(1, T); for t 1:T % 预测从状态转移先验中采样 x f(x) sqrt(Q) * randn(1, N); % 加权用观测似然更新权重 w w .* normpdf(y(t) - g(x), 0, sqrt(R)); w w / sum(w); % 重采样判断 Neff 1 / sum(w.^2); if Neff 0.5 * N idx systematic_resample(w); x x(idx); w ones(1, N) / N; end samples(t, :) x; x_est(t) sum(w .* x); end end function idx systematic_resample(w) % 系统重采样一次随机偏移产生 N 个等距采样点 N numel(w); edges ((0:N-1) rand) / N; edges(end) min(edges(end), 1 - eps); cw cumsum(w); idx zeros(N, 1); j 1; for i 1:N while cw(j) edges(i) j j 1; end idx(i) j; end end使用时先把观测数据准备好再定义转移和观测函数。normpdf来自统计工具箱如果电脑里没装手写exp(-0.5*((y-g(x))/sqrt(R)).^2) / sqrt(2*pi*R)效果完全一样。重采样阈值 0.5*N 是最常见的经验值取 0.3 会更省重采样次数取 0.7 会更激进样本多样性更好但计算量上升。注意每次重采样都会引入额外随机性所以必须保存随机种子这一点在 6.1 节会展开。系统重采样比多项重采样方差更小因为它只抽一个随机数其余 N-1 个点等距排在单位区间上。实现里的edges ((0:N-1) rand) / N就是那个“一次偏移加等距点”用累积分布函数的逆变换把每个区间映射成某个粒子的索引。最后的min(edges(end), 1-eps)是防边界越界最后一个点理论上不可能取到 1但浮点误差下还是兜底一下更稳。3.3 重采样策略与观测噪声很小的退化陷阱重采样选哪种影响的是粒子多样性。多项重采样每次独立抽 N 个序号方差大粒子重复率高系统重采样把 [0,1] 均匀切成 N 段每段取一个点重复率低。实际工程里系统重采样几乎总是更好的默认选项。另一个容易被忽略的是“重采样后要不要加微小抖动”如果需要维持粒子多样性可以给重复粒子叠一层极小幅度的过程噪声但这会引入额外偏差慎用。更隐蔽的坑来自观测噪声 R 的设置。R 设得比真实噪声小一个量级时似然函数变得非常尖权重在第一步就会全部集中到一两个粒子上Neff 瞬间掉到 2 以下。这时候重采样救不回来因为重采样只能复制“还不错的粒子”不能凭空生成“更好的粒子”。现象就是滤波轨迹出现成段的平直粒子全部相同后续预测只在同一条轨迹上加噪声。排查方法很简单把每轮的 Neff 画出来如果在前几步就枯竭先检查 R再检查观测函数是否写对了。状态转移噪声 Q 则相反Q 太小时粒子长期挤在一起多样性不足Q 太大时粒子发散靠似然拉回来需要更多样本。4. 边缘粒子滤波把可解析积分的状态边缘化减少粒子维度4.1 Rao-Blackwellization 为什么能压低方差边缘粒子滤波在英文里常叫 Rao-Blackwellized Particle Filter中文里的“边缘”指的是对部分状态做解析边缘化。做法是把状态拆成两块x_r 保留给粒子x_l 在给定 x_r 和观测的条件下可以解析求解通常是一个卡尔曼滤波问题。粒子只负责采样 x_rx_l 的条件后验用解析高斯分布表达。方差为什么变小可以看条件方差公式Var(X) E[Var(X|Y)] Var(E[X|Y])。用条件期望 E[X|Y] 去替代 X损失的只是条件方差那一项而这部分正是被解析滤波吸收掉了。对比纯粒子滤波RBPF 把一个高维采样问题拆成“低维采样 解析滤波”粒子维度一降Neff 的退化速度大幅放缓。这也是 RBPF 在目标跟踪、同步定位与地图构建里受欢迎的原因位置姿态用粒子地图特征或线性速度部分用卡尔曼各干各的。4.2 粒子采样与卡尔曼更新交替的RBPF主循环假设模型里 x_l 是标量且动态是线性的x_r 通过非线性函数影响 x_l 的演变。为了让结构清楚这里展示主循环的核心代码强调每个粒子都背着一个小卡尔曼滤波器。% RBPF 主循环骨架Np 个粒子P 保存每个粒子对 x_l 的方差 xl_pred zeros(Np, 1); P_pred zeros(Np, 1); for t 1:T % 第 1 步对 x_r 做粒子传播 xr fr(xr) sqrt(Qr) * randn(Np, 1); % 第 2 步对每个粒子做 x_l 的卡尔曼预测 for p 1:Np A a_func(xr(p)); % 线性系数随 x_r 变化 xl_pred(p) A * xl(p); P_pred(p) A^2 * P(p) Ql; end % 第 3 步观测更新与权重更新同时做 for p 1:Np innov y(t) - xl_pred(p); S P_pred(p) R; K P_pred(p) / S; xl(p) xl_pred(p) K * innov; P(p) (1 - K) * P_pred(p); % 权重用边缘似然即观测在高斯预测分布下的密度 w(p) w(p) * normpdf(innov, 0, sqrt(S)); end w w / sum(w); if 1 / sum(w.^2) 0.5 * Np idx systematic_resample(w); xr xr(idx); xl xl(idx); P P(idx); w ones(Np, 1) / Np; end end这段代码的关键在第三层循环中的权重更新权重乘的不是“给定 x_l 后验的似然”而是边际似然即只把 x_r 和观测之间的信息算进来。卡尔曼增益 K 负责修正 x_l而 x_r 的好坏完全由 innovation 的密度体现数值上正好是normpdf(innov, 0, sqrt(S))。如果这一步写成了用完整似然更新权重等于把 x_l 的信息重复计算了一遍滤波会表现得过于自信协方差估计偏小。重采样时最容易漏掉的是 P。xl 被重采样了但 P 还留在旧粒子上下一步的 P_pred 就会和 xl 不匹配卡尔曼增益失真。所以 P 必须跟着 xl 一起按 idx 重排。另一个结构上的经验是能解析的部分尽量拉大。x_r 里如果混进一个弱可观测的维度粒子数需求立刻上升RBPF 的优势就被削弱了。4.3 写RBPF最容易踩的三个代数坑第一状态增广的顺序。RBPF 里 x_r 在粒子中流动但 x_l 的协方差 P 是每个粒子各自维护的矩阵维度必须始终一致。很多人第一步写对了改动模型后 P 的维度忘了跟着改报错位置在卡尔曼增益计算处看起来像“除零”实际上是维度不匹配。建议在进入主循环前用assert(isequal(size(P), [Np, size(xl, 2)]))这类检查把所有数组维度卡死。第二权重更新用边缘似然而不是条件似然。数学上p(y_t | x_{r,1:t}, y_{1:t-1})是高斯分布均值为 xl_pred方差为 S也就是代码里的normpdf(innov, 0, sqrt(S))。如果这里误写成normpdf(y(t) - g([xr(p); xl(p)]), 0, sqrt(R))等于把 x_l 的解析结果也当成了采样变量MCMC 意义上没有错但方差会变大Neff 下降速度明显变快。判断方法是对比两种写法下的 Neff 曲线RBPF 的正确实现曲线应当更平滑。第三可观测性。x_l 的卡尔曼更新依赖观测方程里能看到它。如果某段时间观测对 x_l 完全无信息P 会不降反升权重也退化成只由 x_r 的预测密度决定。这种情况不是代码错而是模型本身在结构上不可观。取舍方法是把不可观的那部分状态挪到 x_r 里或者加一个惩罚先验。5. 调参顺序与选型粒子数、迭代次数与计算成本的取舍5.1 三种方法的适用边界先选结构再选参数对比项变分贝叶斯 (VB)粒子滤波 (PF)边缘粒子滤波 (RBPF)近似对象参数化分布 q(z)后验的经验分布x_r 经验分布 x_l 解析高斯对非线性支持弱依赖变分族灵活性强直接采样只对 x_r 支持非线性典型计算瓶颈迭代轮次 × 状态维度粒子数 × 步长粒子数 × (卡尔曼 积分)初始化敏感度高中中MATLAB 官方函数无通用 VB 函数部分版本工具集有 particleFilter需自行拼装选型逻辑从数据结构出发。如果模型里存在明显可解析的线性高斯块优先考虑 RBPF它的粒子维数最低同样的粒子数下精度最高。如果模型完全非线性且没有高斯结构PF 是兜底方案。VB 适合你要的不是逐时刻的状态轨迹而是参数的完整后验比如高斯混合模型的聚类中心后验它给的是一个分布形状而非样本序列。对做机器学习期末复习或要交课程作业的人来说最简单的判别方式是看关键变量是否连续可导连续且先验共轭用 VB强非线性且转移函数复杂用 PF两者混合用 RBPF。5.2 变分贝叶斯的迭代次数、容差与初值直觉VB 的迭代参数比粒子滤波少但更微妙。先设置一个最大迭代次数比如 300再设置 ELBO 变化量的容差 1e-4 或 1e-6。迭代次数不是越大越好因为均值场分解的偏差不会随迭代消失多跑几千轮只是逼近同一个目标。初值的影响集中在“对称性破缺”的模型上比如高斯混合模型如果初始把两个分量放在同一位置坐标上升可能永远不分开但这在简单的高斯均值模型上没有影响。还有一个从优化工具箱里带过来的习惯用目标函数的变化量做停止判据别直接用参数变化。ELBO 是目标muN 是参数。参数变化很小而 ELBO 还在爬坡说明进入了平坦区域ELBO 先升后降说明某一步期望算错。MATLAB 的优化工具箱在调试时能帮你验证凸问题下的收敛行为但 VB 的更新不是梯度下降不需要也不应该设置学习率。如果看到迭代震荡去检查 q(μ) 更新里是否用了最新的 bN而不是上一轮的旧值。5.3 粒子数与重采样阈值的工程经验值粒子数 N 的经典经验值一维状态 200~1000二维 1000~5000三维以上 5000~20000。这只是一个起点真正要看的指标是重采样频率和 Neff 最低值。如果 Neff 长期低于 0.2N先别急着加粒子而是调整过程噪声 Q 和观测噪声 R。增加 Q 会提高粒子多样性降低 R 会让粒子更集中多数情况下调噪声比重采样阈值有效。在学校实验室搭的机器学习服务器上跑大规模任务时粒子滤波的循环天然可并行每个粒子的预测和加权互不依赖。MATLAB 里可以把第一层 for 换成parfor但要注意重采样步骤必须在并行循环外统一做否则粒子索引会乱。RBPF 的每个粒子有独立的卡尔曼状态同样适合 parfor只要 P、xl、xr 都按切片传进去。阈值 0.5N 和 0.7N 之间的差别远小于噪声参数带来的差别不必在这上面过度纠结先把一套固定的设置跑通再观察 Neff 曲线决定往哪个方向调。6. 在MATLAB里验证这三种实现是否写对的三个技巧6.1 固定随机种子并保存中间状态让每次运行完全可复现粒子滤波和变分贝叶斯都吃随机数。粒子滤波的随机来自初始粒子和重采样VB 如果初始化用了randn同样不可复现。在脚本最上方加一行rng(2024)是最基本的但还不够调试时需要定位到“第几步的重采样出了问题”所以要保存每一步的随机数生成器状态。rng(0); seed_track cell(T, 1); for t 1:T seed_track{t} rng; % 粒子滤波或RBPF的预测与加权步骤 endrng不带参数返回值时取回当前生成器状态记录在 cell 里。下次出错后把生成器恢复成rng(seed_track{t})再单步运行粒子轨迹会完全复现。这个技巧对粒子滤波尤其重要因为重采样一旦混入不同机器或不同 MATLAB 版本产生的随机序列都不同定位到具体哪一步才能把问题从随机性中剥离开。6.2 用ELBO或对数似然曲线定位实现错误VB 实现最常见的错误是更新公式少了一项“‘次要项’曲线不会立即变红而是收敛到错误的 ELBO 值。正确写法是对每轮更新后的 q 分布重新计算完整 ELBO单独存一条曲线。如果 ELBO 最后稳定在某个值但没有全程单调递增说明至少一个因子更新前后没有提高下界问题多半出在该因子更新时用了过期参数。粒子滤波这边可以计算对数似然log p(y_{1:T})的估计公式为每一轮权值和的均值之积取对数。写一个监控变量log_lik log_lik log(mean(unnorm_w))如果它出现异常跳变通常是观测函数写反了或者观测噪声 R 数量级错了。这个方法比看状态轨迹更灵敏因为轨迹误差会被局部平滑掩盖而似然直接把每一步的预测质量摊开。6.3 让工具箱实现当对照互相对比状态估计轨迹MATLAB 较新版本里控制或系统辨识相关工具箱提供particleFilter对象输入which particleFilter可以确认本机是否可用。用它构建同样的模型跑一遍再与自己的 SIR 实现做轨迹对比。两者不应逐点完全一致因为随机性不同但状态估计曲线的趋势和波动幅度应当接近。如果自己的实现输出明显偏离最常见原因是重采样索引方向写错导致粒子被整体倒序重排。另一种对照是用unscentedKalmanFilter做参考线。无迹卡尔曼滤波对弱非线性有一阶精度粒子滤波在粒子数充足时更精确。两者之间的差距应该在合理范围内差距太小说明粒子数可能偏多或问题太线性差距太大优先查重采样和似然函数。最后一招是把过程噪声和观测噪声都设成极小值此时模型近似确定性粒子滤波的输出应该趋近于状态转移函数本身的轨迹。如果这个条件下还发散代码一定存在结构性错误与参数无关。本文还有配套的精品资源点击获取 
