变分贝叶斯、粒子滤波与边缘PF的MATLAB实现与避坑指南

简介这份资源面向机器学习进阶学习者与研究者聚焦变分贝叶斯、粒子滤波及边缘粒子滤波三类概率推断与状态估计方法帮助读者打通从理论推导到MATLAB实现的完整链路。包内共44个文件以27个PDF课件、12个.m脚本为主另含2个pptx讲义、2个ipynb笔记本和1份docx说明压缩包约64.48MB。内容覆盖变分推断、蒙特卡洛、马尔可夫链、隐马尔可夫模型、卡尔曼滤波、高斯混合模型、主题模型与神经网络等主题并配有粒子滤波、重采样、EM、吉布斯采样等可运行代码。已有649人学习下载。读者可借助徐亦达老师课件系统梳理贝叶斯推断脉络通过修改与运行脚本复现算法细节在目标跟踪、动态系统建模等场景中积累调参与排错经验适合作为课程学习或科研入门的配套实践材料。1. 从一组课件和 MATLAB 代码说起变分贝叶斯、粒子滤波与边缘 PF 到底解决什么问题如果你手头正好有一批关于变分贝叶斯、粒子滤波以及边缘粒子滤波的课件和 MATLAB 代码想真正把它们跑起来、用起来而不是停留在翻 PPT 的阶段那这篇笔记就是按这个目标写的。这三个词放在一起不是拼盘变分贝叶斯负责在参数后验上做近似推断粒子滤波负责在非线性非高斯状态空间里做序列蒙特卡洛边缘 PF 则把两者接起来对线性子结构做边缘化、对非线性部分用粒子采样。很多做目标跟踪、导航、时间序列状态估计的工程师最后都会落到这套组合上。MATLAB 在这里的角色很实际——矩阵运算快、绘图直观、调试方便适合把公式一行行验证成可复现的数值结果。下面按“概念立住、代码跑通、参数调对、坑绕开”的顺序展开新手能跟着走熟手能直接看边界条件。2. 变分贝叶斯在 MATLAB 里怎么落地从 ELBO 到可运行脚本2.1 为什么选变分推断而不是直接采样变分贝叶斯Variational Bayes, VB的核心思路是当后验分布 p(Z|X) 没法解析计算时找一个简单的分布 q(Z) 去逼近它把推断问题转成优化问题。优化的目标就是证据下界ELBOELBO E_q[log p(X,Z)] - E_q[log q(Z)]第一项是联合似然的期望第二项是熵。最大化 ELBO 等价于最小化 q 和真实后验之间的 KL 散度。相比 MCMC 采样VB 收敛快、适合在线场景代价是近似有偏q 的族选得不好结果会系统性偏离。在 MATLAB 里做 VB最常见的场景是高斯混合模型GMM的参数推断或者线性高斯状态空间模型的参数学习。选它的理由很直接你不需要等几千次采样收敛迭代几十轮就能看到参数稳定。2.2 用 MATLAB 写一个变分贝叶斯 GMM 的最小实现下面这段代码实现一维数据的变分贝叶斯高斯混合核心是更新每个分量的后验参数均值、精度和混合权重。数据用两个高斯叠加生成方便验证。% 变分贝叶斯高斯混合模型最小实现 % 数据两个一维高斯混合 rng(42); N 500; data [randn(N/2,1)*0.5 2; randn(N/2,1)*0.8 - 3]; K 4; % 初始分量数设大一点让VB自动裁剪 alpha0 1e-3; % 权重Dirichlet先验 beta0 1e-3; % 均值精度先验 nu0 1e-3; % 精度Gamma先验形状 m0 mean(data); % 均值先验中心 W0 1; % 精度先验尺度 % 初始化后验参数 alpha alpha0 * ones(1,K); beta beta0 * ones(1,K); nu nu0 * ones(1,K); m linspace(min(data), max(data), K); W W0 * ones(1,K); maxIter 200; tol 1e-6; ELBO_old -inf; for iter 1:maxIter % E步计算责任度 r_nk E_log_pi psi(alpha) - psi(sum(alpha)); E_log_lambda psi(nu/2) log(2*W) - 1; E_mu_lambda nu .* W beta .* (m.^2); log_rho zeros(N,K); for k 1:K log_rho(:,k) E_log_pi(k) 0.5*E_log_lambda(k) ... - 0.5*(E_mu_lambda(k) - 2*m(k)*data data.^2) * (nu(k)*W(k)); end log_rho log_rho - max(log_rho,[],2); r exp(log_rho) ./ sum(exp(log_rho),2); % M步更新后验参数 Nk sum(r,1) eps; xbar (r * data) ./ Nk; S zeros(1,K); for k 1:K S(k) sum(r(:,k) .* (data - xbar(k)).^2); end alpha alpha0 Nk; beta beta0 Nk; m (beta0*m0 Nk.*xbar) ./ beta; nu nu0 Nk; W 1 ./ (1/W0 S (beta0.*Nk./(beta0Nk)).*(xbar - m0).^2); % 计算ELBO简化版用于监控收敛 ELBO sum(alpha) - K*alpha0 ... 0.5*sum(E_log_lambda .* (nu - 1)) ... sum(Nk .* (E_log_pi - 0.5*E_log_lambda)) ... - 0.5*sum(Nk .* (nu .* W .* (xbar - m).^2 1./beta)); if abs(ELBO - ELBO_old) tol fprintf(收敛于第 %d 轮ELBO %.4f\n, iter, ELBO); break; end ELBO_old ELBO; end % 输出有效分量权重占比 1% weight alpha / sum(alpha); valid weight 0.01; fprintf(有效分量数%d\n, sum(valid)); fprintf(权重%s\n, mat2str(weight(valid), 3)); fprintf(均值%s\n, mat2str(m(valid), 3));逻辑说明E 步用 digamma 函数 psi 计算期望得到每个数据点属于各分量的责任度 rM 步用 r 加权更新 alpha、beta、nu、m、W。ELBO 用来判断收敛不需要精确计算只要单调变化趋势对就行。参数说明K 初始设大VB 的自动相关性判定ARD会把多余分量的权重压到接近零alpha0、beta0、nu0 设很小是给弱先验让数据说话maxIter 和 tol 控制迭代停止。2.3 收敛判断与分量裁剪的实操细节跑上面代码时你会看到 ELBO 在前几十轮快速上升之后变平。如果 ELBO 震荡不收敛通常是两个原因一是 K 设得太大且先验太强分量之间抢数据二是数据尺度没归一化导致数值溢出。我一般会先把数据标准化到零均值单位方差再跑 VB最后把均值和方差还原回原始尺度。分量裁剪看的是权重 alpha/sum(alpha)不是看 Nk。有些分量 Nk 不小但 alpha 被先验拉平说明它和别的分量高度重叠该合并。判断标准权重低于 1% 且均值和其他分量距离小于 0.5 个标准差就可以删。3. 粒子滤波的 MATLAB 实现从重要性采样到重采样3.1 粒子滤波为什么在非线性非高斯下不可替代卡尔曼滤波要求系统线性、噪声高斯。一旦状态转移或观测方程非线性或者噪声是多峰的卡尔曼的假设就崩了。粒子滤波用一组带权重的样本粒子去近似后验分布不依赖线性化也不要求高斯假设。代价是计算量随粒子数线性增长维度高了会遭遇权重退化。标准粒子滤波的流程初始化粒子 → 预测按状态转移采样→ 更新权重按观测似然→ 重采样解决退化→ 估计状态。MATLAB 里实现这套流程关键是向量化操作避免 for 循环逐个粒子更新。3.2 一个可复现的粒子滤波 MATLAB 脚本下面用一维非线性模型做演示状态转移 x_k 0.5x_{k-1} 25x_{k-1}/(1x_{k-1}^2) 8cos(1.2k) 噪声观测 y_k x_k^2/20 噪声。这是经典的 UNGM 模型常用来测试粒子滤波性能。% 一维非线性模型粒子滤波 rng(7); T 50; % 时间步 Np 500; % 粒子数 x zeros(1,T); y zeros(1,T); x(1) 0.1; % 生成真实状态和观测 for k 2:T x(k) 0.5*x(k-1) 25*x(k-1)/(1x(k-1)^2) 8*cos(1.2*k) randn*1; end for k 1:T y(k) x(k)^2/20 randn*1; end % 粒子滤波 particles randn(1,Np)*2; % 初始粒子 weights ones(1,Np)/Np; x_est zeros(1,T); N_eff zeros(1,T); for k 1:T % 预测 particles 0.5*particles 25*particles./(1particles.^2) ... 8*cos(1.2*k) randn(1,Np)*1; % 更新权重 likelihood exp(-0.5*(y(k) - particles.^2/20).^2); weights weights .* likelihood; weights weights / (sum(weights) eps); % 估计 x_est(k) sum(weights .* particles); % 有效粒子数 N_eff(k) 1 / sum(weights.^2); % 重采样系统重采样 if N_eff(k) Np/2 cdf cumsum(weights); cdf(end) 1; u ([0:Np-1] rand)/Np; idx zeros(1,Np); j 1; for i 1:Np while u(i) cdf(j) j j 1; end idx(i) j; end particles particles(idx); weights ones(1,Np)/Np; end end % 绘图对比 figure; plot(1:T, x, k-, LineWidth, 1.5); hold on; plot(1:T, x_est, r--, LineWidth, 1.5); legend(真实状态, PF估计); xlabel(时间步); ylabel(状态值); title(粒子滤波跟踪结果);逻辑说明预测步按状态方程传播粒子更新步用观测似然乘权重并归一化重采样在有效粒子数 N_eff 低于 Np/2 时触发用系统重采样保持多样性。参数说明Np 取 500 是精度和速度的折中一维问题 200 就够四维以上建议 1000 起步重采样阈值 Np/2 是经验值低于这个数权重退化明显过程噪声和观测噪声的方差要匹配实际模型设错会导致粒子发散或过度集中。3.3 重采样策略对比系统重采样、残差重采样与分层重采样系统重采样实现简单但同一权重大的粒子会被复制多次多样性损失快。残差重采样先按权重整数部分复制再对余下部分做多项式重采样多样性稍好。分层重采样把 [0,1] 分成 Np 层每层随机取一点方差比系统重采样小。重采样方法计算复杂度多样性保持适用场景系统重采样O(Np)一般实时性要求高残差重采样O(Np)较好粒子数少分层重采样O(Np)好精度要求高MATLAB 里系统重采样用 cumsum 加二分查找最快残差重采样多一步整数复制分层重采样用 rand 分层。我一般默认用系统重采样粒子数低于 200 时换分层。4. 边缘粒子滤波边缘 PF把线性部分解析消掉4.1 边缘 PF 的核心思想与适用条件边缘粒子滤波Marginalized Particle Filter, MPF针对的是这样一类模型状态可以拆成线性部分和非线性部分线性部分的条件后验能用卡尔曼滤波解析计算。做法是对非线性部分用粒子滤波采样对每个粒子对应的线性部分跑一个卡尔曼滤波最后把线性状态边缘化掉。这样做的好处是降维。原本需要在高维空间撒粒子现在只在非线性维度撒线性维度用卡尔曼解析处理粒子数需求大幅下降。适用条件状态方程和观测方程对线性子状态是线性的噪声高斯非线性子状态可以是非线性非高斯。4.2 边缘 PF 的 MATLAB 实现框架下面用一个简化模型演示状态 [x_linear; x_nonlinear]线性部分做匀速运动非线性部分做非线性变换。观测同时包含两部分。% 边缘粒子滤波简化实现 rng(11); T 40; Np 200; % 真实状态 x_lin zeros(2,T); % [位置; 速度] x_non zeros(1,T); x_lin(:,1) [0; 1]; x_non(1) 0.5; for k 2:T x_lin(:,k) [1 1; 0 1] * x_lin(:,k-1) [0.5; 1]*randn*0.1; x_non(k) 0.6*x_non(k-1) 2*x_non(k-1)/(1x_non(k-1)^2) randn*0.3; end y zeros(1,T); for k 1:T y(k) x_lin(1,k) x_non(k)^2/10 randn*0.5; end % 边缘PF particles_non randn(1,Np)*0.5; weights ones(1,Np)/Np; % 每个粒子配一个卡尔曼滤波器 kf_mean repmat([0;1], 1, Np); kf_cov repmat(eye(2)*0.1, [1,1,Np]); x_est zeros(1,T); x_non_est zeros(1,T); A [1 1; 0 1]; Q [0.5;1]*[0.5 1]*0.01; H [1 0]; for k 1:T % 非线性粒子预测 particles_non 0.6*particles_non 2*particles_non./(1particles_non.^2) randn(1,Np)*0.3; % 每个粒子的卡尔曼预测 for i 1:Np kf_mean(:,i) A * kf_mean(:,i); kf_cov(:,:,i) A * kf_cov(:,:,i) * A Q; end % 权重更新用观测残差 for i 1:Np y_pred H * kf_mean(:,i) particles_non(i)^2/10; S H * kf_cov(:,:,i) * H 0.25; K kf_cov(:,:,i) * H / S; innov y(k) - y_pred; kf_mean(:,i) kf_mean(:,i) K * innov; kf_cov(:,:,i) (eye(2) - K*H) * kf_cov(:,:,i); weights(i) weights(i) * exp(-0.5 * innov^2 / S) / sqrt(S); end weights weights / (sum(weights) eps); % 估计 x_est(k) sum(weights .* kf_mean(1,:)); x_non_est(k) sum(weights .* particles_non); % 重采样 if 1/sum(weights.^2) Np/2 cdf cumsum(weights); cdf(end) 1; u ([0:Np-1] rand)/Np; idx zeros(1,Np); j 1; for i 1:Np while u(i) cdf(j), j j 1; end idx(i) j; end particles_non particles_non(idx); kf_mean kf_mean(:,idx); kf_cov kf_cov(:,:,idx); weights ones(1,Np)/Np; end end figure; subplot(2,1,1); plot(1:T, x_lin(1,:), k-, 1:T, x_est, r--); legend(真实位置, MPF估计); title(线性部分); subplot(2,1,2); plot(1:T, x_non, k-, 1:T, x_non_est, r--); legend(真实非线性状态, MPF估计); title(非线性部分);逻辑说明非线性粒子按非线性方程预测每个粒子维护一个卡尔曼滤波器对线性部分做预测和更新权重用观测残差的似然更新重采样时粒子、卡尔曼均值和协方差同步重采样。参数说明Q 是线性部分过程噪声协方差设太小会导致卡尔曼增益过低、跟踪滞后Np 在边缘 PF 里可以比标准 PF 少因为线性维度被解析处理了重采样阈值同样用 Np/2。4.3 边缘 PF 与标准 PF 的计算量对比维度标准 PF 粒子数边缘 PF 粒子数计算量比2 维500200约 2.5:14 维2000300约 6:16 维10000500约 20:1维度越高边缘 PF 的优势越明显。但前提是线性子结构确实存在强行拆一个非线性部分当线性处理卡尔曼滤波会发散整个边缘 PF 就崩了。5. 避坑与排查课件代码跑不通时先看这几条5.1 现象粒子滤波估计值发散和真实状态越差越远原因过程噪声方差设得远小于实际粒子云收缩到一点观测似然再大也拉不回来。或者重采样太频繁粒子多样性几轮就没了。解决先检查过程噪声是否匹配状态转移的实际扰动。用 N_eff 曲线判断如果 N_eff 在重采样后立刻掉到 Np/2 以下说明重采样阈值太松改成 Np/3 或 Np/4。另外重采样后给粒子加一点抖动jitter方差取过程噪声的 1/10。5.2 现象变分贝叶斯 ELBO 不收敛来回震荡原因数据没归一化不同维度尺度差几个数量级ELBO 里各项量级不匹配。或者 K 设得太大分量之间反复抢数据。解决先做 z-score 标准化跑完 VB 再还原。K 从 2 开始试逐步加看 ELBO 是否稳定上升。如果加 K 后 ELBO 反而下降说明先验太强把 alpha0 调小一个数量级。5.3 现象边缘 PF 里卡尔曼滤波协方差矩阵非正定报错原因数值误差累积P (I-KH)P 在多次迭代后失去对称正定性。或者 Q 设得太小P 收缩到接近奇异。解决每步更新后做对称化 P (PP)/2并加一个小的对角正则项 1e-8*eye(2)。Q 不要设为零矩阵至少给 1e-6 量级的对角。5.4 现象MATLAB 跑粒子滤波时内存爆掉原因粒子数设太大或者把每个粒子的历史轨迹都存下来了。T1000、Np10000 时一个 T×Np 的矩阵就是 80MB多个中间变量叠加就爆。解决只存当前时刻的粒子历史轨迹用估计值存。粒子数按维度定一维 200、二维 500、四维 1000超过四维考虑边缘 PF 或 UKF。5.5 现象课件里的代码复制过来报未定义函数原因课件代码依赖作者自己写的工具函数比如系统重采样、digamma 的替代实现课件里没附全。解决MATLAB 自带 psi 函数就是 digamma直接用。重采样函数自己写上面第 3 章的系统重采样代码可以直接抄。遇到其他未定义函数先搜 MATLAB 官方文档有没有同名内置函数没有再自己实现。6. 进阶技巧用边缘 PF 做在线参数学习与模型验证把边缘 PF 从状态估计扩展到参数学习思路是把静态参数增广到非线性状态里用粒子滤波同时估计状态和参数。但参数不随时间变化直接增广会导致粒子退化——因为参数维度的粒子没有过程噪声几轮重采样后全变成同一个值。常见做法是给参数粒子加一个人工抖动抖动方差随迭代次数衰减这就是 Liu-West 重采样。具体操作在边缘 PF 的重采样步骤后对参数粒子做一次收缩变换。设参数粒子为 θ_i均值为 θ_bar收缩公式 θ_i sqrt(1-h^2)(θ_i - θ_bar) θ_bar hrandn*σ_θ其中 h 取 0.01 到 0.1σ_θ 是参数先验标准差。这样既保持了参数粒子的多样性又不会让参数漂移太远。验证方法用仿真数据跑一遍看参数估计是否收敛到真值附近。如果参数估计在真值附近震荡但不收敛说明 h 太大如果参数估计几乎不动说明 h 太小。我一般从 h0.05 开始试看参数轨迹的收敛速度。另一个技巧是用边缘 PF 的似然输出做模型选择。每个模型的边缘似然可以用粒子权重的乘积累加估计比较不同模型的边缘似然选大的那个。这个在 MATLAB 里实现就是多跑几组边缘 PF每组记录权重和最后比较总和。最后说一个我踩过的坑边缘 PF 里卡尔曼滤波的观测矩阵 H 必须和粒子部分的观测方程对齐如果观测方程里线性部分和非线性部分有交叉项比如 y x_lin * x_non那就不能直接拆得先把交叉项归到非线性部分或者对观测方程做线性化。这个对齐错了卡尔曼滤波的残差计算就是错的整个权重更新全废。我一般会在代码里单独写一个观测函数输入线性状态和非线性状态输出预测观测值确保预测和真实观测的维度、量纲一致。希望帮到你。本文还有配套的精品资源点击获取