简介这份资源面向机器学习进阶学习者与研究者聚焦变分贝叶斯、粒子滤波及边缘粒子滤波三类概率推断与状态估计方法帮助读者打通从理论推导到MATLAB实现的关键环节。包内共44个文件以27个PDF课件、12个MATLAB脚本为主另含2个PPTX讲义、2个Jupyter Notebook和1份说明文档压缩包约64.48MB覆盖变分推断、蒙特卡洛、隐马尔可夫模型、卡尔曼滤波、主题模型与神经网络等主题。资源包含徐亦达老师机器学习课件理论讲解与配套代码相互印证可运行vb_normal_gamma、particle_filter、gibbs、kalman_demo等脚本复现算法流程并借助notebook理解中文餐馆过程等非参数模型。目前已有649人学习下载适合希望系统掌握贝叶斯推断与动态系统估计、并积累可复用MATLAB代码的读者。1. 从一组课件和 MATLAB 代码说起变分贝叶斯、粒子滤波与边缘 PF 到底解决什么问题如果你手头正好拿到一份标题写着「机器学习变分贝叶斯、粒子滤波及边缘 PF内容包括大量课件、MATLAB 代码」的资料包第一反应大概率是东西挺全但不知道从哪下手。变分贝叶斯、粒子滤波、边缘粒子滤波这三个词放在一起其实指向的是同一类问题——状态估计与参数推断只不过一个走的是确定性近似路线一个走的是蒙特卡洛采样路线边缘 PF 则是把两者缝在一起的混合方案。它们最典型的落地场景是目标跟踪、导航定位、时间序列预测、通信信道估计这类「系统有隐状态、观测带噪声」的活儿。MATLAB 在这里的角色很实在矩阵运算快、绘图直观、调试方便特别适合把公式一行行翻译成能跑的实验。这篇笔记就按「概念立住 → 最小复现 → 参数怎么调 → 坑在哪」的顺序把这三个东西讲成能照着做的方案适合刚接触状态估计的新手也适合想从纯理论转向工程实现的熟手。2. 变分贝叶斯把难解的后验变成可优化的下界2.1 为什么需要变分推断而不是直接算后验贝叶斯推断的核心目标是求后验分布 $p(z|x)$其中 $z$ 是隐变量$x$ 是观测。理论上用贝叶斯公式就能写出来但分母那个归一化积分 $p(x)\int p(x|z)p(z)dz$ 在高维空间里几乎算不动。传统做法是 MCMC 采样精度高但慢每来一批新数据都要重新跑。变分贝叶斯换了个思路不直接求后验而是找一个简单的分布 $q(z)$ 去逼近它把「推断」变成「优化」。衡量逼近好坏用的是 KL 散度 $KL(q(z)|p(z|x))$但这个量本身也依赖后验没法直接算。于是做一步推导得到证据下界 ELBO$$\mathcal{L}(q)\mathbb{E}{q}[\log p(x,z)]-\mathbb{E}{q}[\log q(z)]$$最大化 ELBO 等价于最小化 KL 散度而且 ELBO 只涉及联合分布 $p(x,z)$ 和变分分布 $q(z)$都是可算的。这就是变分贝叶斯能落地的根本原因——把不可解的积分问题转成了可解的优化问题。选 $q(z)$ 的时候最常用的是平均场假设即把 $q(z)$ 拆成各维独立的乘积 $q(z)\prod_i q_i(z_i)$。这个假设牺牲了变量间的相关性但换来的是坐标上升式的迭代更新每一维都有闭式解。常见做法是共轭先验配共轭似然比如高斯配高斯更新公式直接是均值和精度的代数运算。2.2 用 MATLAB 跑通一个变分贝叶斯高斯混合的最小例子下面这段代码实现的是变分贝叶斯高斯混合VB-GMM的核心迭代用一维数据演示方便看清每一步在更新什么。数据自己造两个高斯簇。% VB-GMM 最小复现一维数据两个簇 rng(1); data [randn(200,1)*0.5 2; randn(200,1)*0.8 - 3]; % 两簇数据 N length(data); K 2; % 假设簇数VB 会自动压制多余簇 % 先验参数共轭高斯-逆Gamma mu0 0; beta0 1e-3; a0 1e-3; b0 1e-3; % 初始化变分参数 mu linspace(min(data), max(data), K); beta ones(K,1); a ones(K,1); b ones(K,1); alpha ones(K,1); % Dirichlet 参数 maxIter 200; tol 1e-6; Lold -inf; for iter 1:maxIter % E 步计算责任度 r_nk log_rho zeros(N,K); for k 1:K E_log_lambda psi(a(k)) - log(b(k)); E_log_pi psi(alpha(k)) - psi(sum(alpha)); log_rho(:,k) E_log_pi 0.5*E_log_lambda ... - 0.5*log(2*pi) - 0.5*(a(k)/b(k)) ... .* ((data-mu(k)).^2 1/beta(k)); end log_rho log_rho - max(log_rho,[],2); % 数值稳定 r exp(log_rho) ./ sum(exp(log_rho),2); Nk sum(r,1); % M 步更新变分参数 xbar (r*data) ./ Nk; for k 1:K beta(k) beta0 Nk(k); mu(k) (beta0*mu0 Nk(k)*xbar(k)) / beta(k); a(k) a0 Nk(k)/2; b(k) b0 0.5*sum(r(:,k).*(data-mu(k)).^2) ... 0.5*beta0*Nk(k)/beta(k)*(xbar(k)-mu0)^2; alpha(k) 1 Nk(k); end % 计算 ELBO 判断收敛 L 0; for k 1:K L L gammaln(a0) - a0*log(b0) (a0-1)*(psi(a(k))-log(b(k))) ... - b0*a(k)/b(k) 0.5*log(beta0) - 0.5*log(beta(k)) ... - 0.5*beta0/beta(k) - 0.5*beta0*(mu(k)-mu0)^2 ... gammaln(a(k)) - (a(k)-1)*psi(a(k)) a(k) ... gammaln(1) - gammaln(alpha(k)); end L L sum(sum(r .* log(r eps))); if abs(L - Lold) tol, break; end Lold L; end fprintf(迭代 %d 次收敛簇均值: %.3f, %.3f\n, iter, mu(1), mu(2));这段代码的逻辑分三块。E 步算每个数据点属于每个簇的责任度 $r_{nk}$用的是变分参数下的期望对数似然注意psi是 digamma 函数对应 Gamma 分布期望的对数。M 步按共轭更新公式刷新mu、beta、a、b、alpha其中beta是精度a/b控制方差。最后算 ELBO 判断收敛log_rho减最大值是为了防止指数溢出这是数值计算里的常规操作。参数上要留意几个点。beta0和a0、b0设得很小1e-3是弱先验让数据说话如果先验太强簇会被强行拉到初始值附近。alpha初始化为 1VB 的一个好处是多余簇的Nk会趋近 0alpha随之收缩自动实现模型选择不用像 EM 那样手动定 K。maxIter给 200 足够一般几十次就收敛。如果发现两个簇均值几乎重合多半是初始化太近把mu的初始值拉开一点。2.3 变分贝叶斯的适用边界与常见误用变分贝叶斯不是万能的。平均场假设丢掉了后验相关性方差通常被低估做不确定性量化时要心里有数。如果后验是多峰的平均场容易卡在某个峰附近这时候要么换更灵活的变分族要么退回 MCMC。另一个常见误用是把 ELBO 当成绝对指标去比较不同模型ELBO 只在同一模型内比较才有意义跨模型比大小没有理论保证。在工程里变分贝叶斯最值得投入的场景是在线学习和大规模数据每来一批数据做一次坐标上升比 MCMC 重新采样快得多。如果数据量小、对精度要求极高MCMC 仍然是更稳的选择。3. 粒子滤波用带权样本逼近任意后验3.1 从重要性采样到序贯重要性重采样粒子滤波解决的是非线性、非高斯状态空间模型里的状态估计问题。模型写成$$z_t f(z_{t-1}) u_t,\quad x_t h(z_t) v_t$$其中 $f$ 和 $h$ 都可以是非线性的$u_t$、$v_t$ 可以是非高斯噪声。卡尔曼滤波在这种场景下直接失效因为高斯假设不成立。粒子滤波的思路是用一组带权粒子 ${z_t^{(i)}, w_t^{(i)}}_{i1}^N$ 去近似后验分布粒子密集的地方概率高。核心是序贯重要性采样从建议分布 $q(z_t|z_{t-1},x_t)$ 里采样按重要性权重 $w_t \propto w_{t-1} \cdot p(x_t|z_t)p(z_t|z_{t-1})/q(z_t|z_{t-1},x_t)$ 更新。最常用的建议分布就是状态转移先验 $p(z_t|z_{t-1})$此时权重简化为似然 $p(x_t|z_t)$实现最简单。但纯重要性采样跑几轮后权重会退化几乎所有粒子权重趋近 0只剩一个粒子权重为 1。解决办法是重采样按权重重抽粒子把权重拉回均匀。重采样会带来样本贫化所以通常加一个有效粒子数 $N_{eff}1/\sum (w^{(i)})^2$ 的判断低于阈值才重采样。3.2 MATLAB 实现一个可复现的粒子滤波跟踪器下面用一维非线性模型演示状态转移 $z_t 0.5z_{t-1} 25z_{t-1}/(1z_{t-1}^2) 8\cos(1.2t)$观测 $x_t z_t^2/20$这是粒子滤波文献里的经典测试模型。% 粒子滤波最小复现一维非线性模型 rng(2); T 50; N 500; % 时间步、粒子数 z zeros(T,1); x zeros(T,1); z(1) 0.1; for t 2:T z(t) 0.5*z(t-1) 25*z(t-1)/(1z(t-1)^2) 8*cos(1.2*t) randn; end x z.^2/20 randn(T,1)*0.5; % 观测噪声 % 粒子滤波 particles randn(N,1)*2; % 初始粒子 weights ones(N,1)/N; z_est zeros(T,1); z_est(1) mean(particles); for t 2:T % 预测按状态转移采样 particles 0.5*particles 25*particles./(1particles.^2) ... 8*cos(1.2*t) randn(N,1); % 更新按似然算权重 likelihood exp(-0.5*((x(t) - particles.^2/20).^2)/0.5^2); weights weights .* likelihood; weights weights / sum(weights); % 估计 z_est(t) sum(weights .* particles); % 重采样判断 Neff 1 / sum(weights.^2); if Neff N/2 idx randsample(N, N, true, weights); particles particles(idx); weights ones(N,1)/N; end end rmse sqrt(mean((z_est - z).^2)); fprintf(粒子滤波 RMSE: %.3f\n, rmse);逻辑上分预测、更新、重采样三步。预测步按状态转移方程把粒子往前推加过程噪声。更新步用观测似然给每个粒子打分权重归一化。估计值取加权平均。重采样用randsample按权重重抽Neff N/2是常用阈值。参数上N是最关键的粒子数太少估计抖动大太多计算慢。一维问题 500 够用高维问题要上千甚至上万。观测噪声标准差 0.5 要和实际匹配设小了权重会过度集中设大了粒子区分度不够。重采样阈值N/2是经验值也有用2N/3的阈值越低重采样越频繁样本贫化风险越高。3.3 粒子数、建议分布与重采样策略怎么选粒子数选择有个粗略经验状态维度每增加一维粒子数大致要翻倍才能维持同样的估计精度这就是维数灾难。高维场景下标准粒子滤波基本不可用需要换边缘粒子滤波或者辅助粒子滤波。建议分布的选择直接影响效率。用状态转移先验最简单但如果观测信息量大大量粒子会落在低似然区浪费算力。改进办法是用观测信息构造更好的建议分布比如扩展卡尔曼滤波或无迹卡尔曼滤波生成建议分布这就是 EKF-PF 和 UKF-PF。代价是要算雅可比矩阵或 sigma 点实现复杂度上升。重采样策略上多项式重采样最简单但方差大系统重采样和分层重采样方差更小残差重采样在并行实现里友好。MATLAB 里randsample默认是多项式重采样对精度要求高的场景可以自己写系统重采样。4. 边缘粒子滤波把变分贝叶斯和粒子滤波缝在一起4.1 边缘 PF 要解决的核心矛盾标准粒子滤波在高维状态空间里会崩因为粒子要覆盖整个空间。但很多实际问题里状态可以拆成两部分一部分是线性高斯的可以用卡尔曼滤波精确处理另一部分是非线性非高斯的才需要粒子。边缘粒子滤波Marginalized Particle Filter也叫 Rao-Blackwellized PF就是利用这个结构把线性部分解析地边缘化掉只对非线性部分采样。数学上假设状态 $z_t [z_t^l, z_t^n]$其中 $z_t^l$ 在给定 $z_t^n$ 时条件线性高斯。那么后验可以分解$$p(z_{0:t}^l, z_{0:t}^n|x_{1:t}) p(z_{0:t}^l|z_{0:t}^n,x_{1:t}) \cdot p(z_{0:t}^n|x_{1:t})$$第一项用卡尔曼滤波解析计算第二项用粒子滤波。每个粒子携带一个卡尔曼滤波器的均值和协方差。这样粒子只需要覆盖非线性部分的低维空间效率大幅提升。边缘 PF 和变分贝叶斯的联系在于当线性部分的参数也未知时可以用变分贝叶斯在线估计这些参数再喂给卡尔曼滤波。这就是标题里三个东西放在一起的工程意义——变分贝叶斯管参数卡尔曼管线性状态粒子管非线性状态。4.2 一个可运行的一维边缘 PF 骨架下面用简化模型演示边缘 PF 的结构状态 $z_t [z_t^l, z_t^n]$$z_t^l$ 线性$z_t^n$ 非线性。% 边缘 PF 骨架线性部分卡尔曼非线性部分粒子 rng(3); T 40; N 200; % 真实状态 zl zeros(T,1); zn zeros(T,1); zl(1) 0; zn(1) 0.5; for t 2:T zl(t) 0.9*zl(t-1) 0.1*zn(t-1) randn*0.3; zn(t) 0.5*zn(t-1) 25*zn(t-1)/(1zn(t-1)^2) randn*0.5; end x zl zn.^2/10 randn(T,1)*0.4; % 观测 % 边缘 PF particles_n randn(N,1)*0.5; % 非线性粒子 weights ones(N,1)/N; kalman_mu zeros(N,1); kalman_P ones(N,1); % 每个粒子的卡尔曼状态 Q 0.3^2; R 0.4^2; z_est zeros(T,1); for t 2:T for i 1:N % 非线性粒子预测 particles_n(i) 0.5*particles_n(i) 25*particles_n(i)/(1particles_n(i)^2) randn*0.5; % 卡尔曼预测线性部分依赖非线性粒子 kalman_mu(i) 0.9*kalman_mu(i) 0.1*particles_n(i); kalman_P(i) 0.9^2*kalman_P(i) Q; % 卡尔曼更新 K kalman_P(i) / (kalman_P(i) R); kalman_mu(i) kalman_mu(i) K*(x(t) - kalman_mu(i) - particles_n(i)^2/10); kalman_P(i) (1-K)*kalman_P(i); % 权重更新 weights(i) weights(i) * exp(-0.5*(x(t)-kalman_mu(i)-particles_n(i)^2/10)^2/R); end weights weights / sum(weights); z_est(t) sum(weights .* (kalman_mu particles_n.^2/10)); % 重采样 if 1/sum(weights.^2) N/2 idx randsample(N, N, true, weights); particles_n particles_n(idx); kalman_mu kalman_mu(idx); kalman_P kalman_P(idx); weights ones(N,1)/N; end end fprintf(边缘 PF 完成末步估计: %.3f\n, z_est(end));这段代码的关键在于每个粒子内部跑一个卡尔曼滤波器。非线性粒子负责预测和采样卡尔曼负责在线性部分做最优估计。权重更新用的是卡尔曼预测的观测似然。重采样时粒子的卡尔曼状态要一起重抽否则会错位。参数上Q和R是卡尔曼的命脉设错了整个估计会偏。Q反映线性状态的过程噪声R反映观测噪声一般用离线数据估或者自适应调。粒子数N比标准 PF 可以少很多因为非线性维度低。如果非线性部分也高维边缘 PF 的优势会减弱。4.3 变分贝叶斯在边缘 PF 里的接入点当线性部分的噪声协方差 $Q$、$R$ 未知时标准做法是人工调参费时费力。变分贝叶斯可以在这里接入把 $Q$、$R$ 当成隐变量给逆 Gamma 先验在每个时间步用变分更新估计它们。这样卡尔曼滤波用的是在线估计的噪声参数适应非平稳环境。接入方式是在每个粒子的卡尔曼更新后加一步变分参数更新。更新公式和 2.2 节类似只是数据变成卡尔曼的预测误差。这一步计算量不大但能显著提升鲁棒性。要注意的是变分更新和粒子权重的耦合如果变分收敛太慢权重会滞后实践中一般给变分更新设一个学习率或者限制每步更新幅度。5. 避坑与排查这几个地方最容易翻车5.1 粒子权重全变成 NaN现象跑几十步后权重全是 NaN估计值消失。原因通常是似然计算里出现了exp下溢或者0/0。观测噪声设得太小似然值极小exp直接变 0归一化时除以 0。解决在算似然前减最大值或者用对数域计算权重再归一化。MATLAB 里可以写log_w log_w - max(log_w)再exp。5.2 变分贝叶斯簇数自动收缩到 1现象VB-GMM 跑完只剩一个簇明明数据有两簇。原因多半是初始化时mu太接近或者先验alpha设得太小导致 Dirichlet 过度稀疏。解决把初始mu拉开alpha初始化给大一点比如 1 而不是 1e-3或者多跑几次不同初始化取 ELBO 最高的。5.3 边缘 PF 估计值系统性偏移现象估计曲线整体偏离真实值不是随机抖动而是固定偏差。原因通常是卡尔曼部分的观测方程写错了或者非线性粒子和线性状态的耦合项漏了。边缘 PF 里观测往往同时依赖两部分状态更新时要把耦合项算进去。解决拿一步数据手算一遍卡尔曼更新和代码输出对确认观测矩阵和耦合项没错。5.4 重采样后估计精度反而下降现象加了重采样RMSE 比不加还大。原因是重采样太频繁导致样本贫化粒子多样性丢失。解决调低重采样阈值或者用残差重采样、系统重采样替代多项式重采样。另一个办法是重采样后加一点抖动正则化粒子滤波给粒子加小噪声恢复多样性。5.5 MATLAB 版本兼容导致的函数报错现象换台机器跑randsample或psi报未定义。原因randsample属于 Statistics and Machine Learning Toolboxpsi属于核心但旧版本行为有差异。解决确认工具箱装了没有ver命令查一下。没有工具箱就自己写重采样psi可以用digamma替代或者查表。课件里的代码如果跑不通先查依赖函数别急着改算法。6. 把三个方法串成一条验证链从单步手算到端到端对比学这三个东西最容易犯的错是「代码跑通了但不知道对不对」。我的习惯是搭一条验证链从最小可验证单元往上走。第一步变分贝叶斯先拿一维高斯数据手算一遍 ELBO确认更新公式和代码一致。第二步粒子滤波拿线性高斯模型跑结果应该和卡尔曼滤波接近差太多说明权重或重采样有问题。第三步边缘 PF 拿线性模型跑此时非线性部分退化成线性边缘 PF 应该和标准卡尔曼几乎重合这是最强的正确性检验。具体操作上可以写一个统一的测试脚本用同一组数据分别跑卡尔曼、标准 PF、边缘 PF对比 RMSE 和运行时间。下面这个表格是我常用的对比维度方法适用模型一维 RMSE 参考计算耗时相对值调参难度卡尔曼滤波线性高斯最低1低标准粒子滤波非线性非高斯中50中边缘粒子滤波部分线性部分非线性接近卡尔曼10高变分贝叶斯参数未知依赖近似质量5中验证的时候有个技巧固定随机种子让三个方法看到完全相同的观测序列这样对比才公平。MATLAB 里rng(seed)放在数据生成之前每个方法跑之前重置一次保证噪声序列一致。另一个我踩过的坑是拿边缘 PF 去跑纯非线性模型结果还不如标准 PF。边缘 PF 的优势建立在「存在条件线性子结构」这个前提上如果模型没有这个结构强行拆只会增加复杂度。判断方法很简单看状态转移和观测方程能不能写成对某部分状态线性的形式不能就别用。最后说一个调试习惯。粒子滤波和边缘 PF 出问题时别盯着最终估计曲线看先把中间量打出来权重分布、有效粒子数、卡尔曼协方差。权重直方图如果大部分接近 0说明建议分布太差有效粒子数一直低于阈值说明重采样太频繁卡尔曼协方差发散说明 Q 或 R 设错了。这些中间量比最终结果更能定位问题。希望帮到你。本文还有配套的精品资源点击获取
