变分贝叶斯、粒子滤波与边缘粒子滤波：从原理到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) 还是绕不开后验的归一化常数，于是把目标等价改写成最大化证据下界（ELBO）：

ELBO(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，想让先验弱就同时取小值。初始值给 muN=mean(y)、lamN=1/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_r，x_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 三种方法的适用边界：先选结构再选参数

选型逻辑从数据结构出发。如果模型里存在明显可解析的线性高斯块，优先考虑 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的预测与加权步骤 end

rng不带参数返回值时，取回当前生成器状态，记录在 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做参考线。无迹卡尔曼滤波对弱非线性有一阶精度，粒子滤波在粒子数充足时更精确。两者之间的差距应该在合理范围内：差距太小说明粒子数可能偏多或问题太线性；差距太大优先查重采样和似然函数。最后一招是把过程噪声和观测噪声都设成极小值，此时模型近似确定性，粒子滤波的输出应该趋近于状态转移函数本身的轨迹。如果这个条件下还发散，代码一定存在结构性错误，与参数无关。

本文还有配套的精品资源，点击获取