Copula变分贝叶斯与信息几何：突破均值场假设的MATLAB实现

简介这是一份面向机器学习和统计推断研究者的 Matlab 项目源码完整复现了论文“Copula Variational Bayes inference via information geometry”中的算法。其核心是用 Copula 描述随机变量间复杂的非线性依赖再借助信息几何在概率分布空间上优化变分贝叶斯近似特别适合理解贝叶斯网络、混合模型中的近似推断与不确定性量化也适合有一定贝叶斯统计背景、希望动手实现现代变分推断算法的读者。压缩包共 42 个文件、约 2.91MB以 25 个 .m 脚本为主体覆盖双变量高斯、高斯混合、蒙特卡洛对比等实验模块另有 .md 说明文档、.png/.eps 示意图和 .html 预览页目录还包含 README 与主运行脚本便于按流程复现。当前已有 242 人学习。代码中能直接看到 ELBO 推导、Copula 参数化、变分更新、EM 聚类评估及绘图等关键函数既能辅助论文复现也可替换数据或修改参数后用于自己的贝叶斯建模实验。1. 变分推断里最该被怀疑的假设均值场独立性正是 CVB 要拆掉的墙变分贝叶斯Variational Bayes最常用的近似路径是 mean-field把后验分布拆成若干个独立因子的乘积再用坐标上升法迭代求解 ELBO。计算是快但代价很直接——变量之间的相关性被整个丢掉了后验分布越“歪”近似结果越不像。Copula Variational Bayes 换了一条路边缘分布仍然用简单分布近似变量间的依赖结构单独交给 copula 密度函数建模而论文题目里的 information geometry 则给出了如何在概率分布空间中度量这个近似的好坏。这份 MATLAB 工程就是论文“Copula Variational Bayes inference via information geometry”算法的一份完整落地实现包含双高斯和混合高斯两套从参数设置、迭代求解到绘图验证的流程适合正在调试变分推断相关问题、想绕开均值场假设去建模后验相关性的研究生和工程师。2. 先跑起来从目录结构到 A0_MAIN 双高斯示例的完整复现路径2.1 目录结构两个独立场景两套 README别跑错目录解压之后不要急着双击任何一个 .m 文件先把这个代码包的目录结构看清楚。整个工程按论文里的两类实验场景分成两大块CVB for bivariate Gauss双变量高斯场景和CVB for Gaussian mixture高斯混合场景。两个场景各自带一个 README.md前者的主入口是A0_MAIN_biGauss.m后者的主入口是A0_MAIN.m文件名风格完全不同一眼就能分辨。容易踩的一个坑是两个场景共用了一部分名称相同的函数比如Func_contourGauss.m和gridxy.m在两个目录里都存在。如果图省事用addpath把两个目录一次性全部加进 MATLAB 路径运行时会因为同名文件覆盖而出现诡异的行为。我的习惯是把当前目录用cd切到目标场景的子文件夹里不依赖路径搜索这样最保险。下面这张表列出的是两个场景里承担关键职责的文件方便对照着看文件在流程里的角色所属场景A0_MAIN_biGauss.m双高斯场景总入口负责参数初始化、调用迭代、触发绘图bivariate GaussB0_CVB.mCVB 迭代主体承载变分参数的循环更新逻辑bivariate GaussFunc_CVB_biGauss.m双高斯模型的 CVB 单步更新函数核心计算在这里bivariate Gausssigma2SIGMA.m把边缘方差向量 σ² 组装成完整协方差矩阵 Σbivariate GaussC0_PLOT.m生成真实后验与 CVB 近似的等高线对比图bivariate GaussFunc_KLDMultiGauss.m计算两个多元高斯分布之间的 KL 散度两场景共用gridxy.m在指定坐标范围绘制均匀网格配合等高线显示两场景共用A0_MAIN.m混合高斯场景总入口Gaussian mixtureB0_SETTING.m混合分量数、先验超参、迭代次数等全部集中在这里Gaussian mixtureC0_RUN.m调用混合场景的主迭代循环Gaussian mixtureFunc_CVB123.m混合高斯场景的 CVB 核心更新函数Gaussian mixtureFunc_VB.m标准 mean-field VB 对照实现用于和 CVB 做效果对比Gaussian mixturec2_MonteCarlo.m蒙特卡洛重复实验的流水线脚本Gaussian mixturec1_para_groundtruth.m定义 ground truth 参数供数据生成和评估使用Gaussian mixture另外根目录下还有一个intro文件夹里面有Fig2_pythagore.png、VAE_vs_CVB_0.png和index.html对应论文里的概念图和 VAE 与 CVB 的对比示意。读代码之前建议先打开这两张图它们对理解信息几何视角下 copula 的作用位置非常有帮助。2.2 复现第一步跑通双高斯示例我用的是 MATLAB R2023b这个工程只用了基础数值计算功能不依赖任何额外工具箱旧版本同样能跑。打开 MATLAB把当前目录切到CVB for bivariate Gauss然后直接执行主入口脚本% 切换到双高斯场景目录 cd(你的本地路径/Copula-Variational-Bayes-master/CVB for bivariate Gauss); % 运行主入口 A0_MAIN_biGauss这个脚本做的事一句话概括就是生成或读入一个已知参数的双变量高斯后验用 CVB 去近似它最后用C0_PLOT.m画一张对比图。运行结束后Figures目录下应该会出现一张等高线图包含真实后验、CVB 近似和标准 VB 近似三条等高线。如果不想直接跑黑盒主函数可以手动一步步执行来理解每个环节的输入输出。下面这段代码展示了双高斯场景最核心的操作顺序% 手动复现双高斯 CVB 的关键步骤 mu_true [0 0]; % 真实后验均值二维 Sigma_true [1 0.6; 0.6 1]; % 真实后验协方差相关系数 0.6 q_mu [0.1 -0.2]; % 变分分布初始均值 q_sigma2 [1.2 1.0]; % 变分分布初始边缘方差 % 用辅助函数把边缘方差拼成协方差矩阵 Sigma_q sigma2SIGMA(q_sigma2); % 单步更新输入当前变分参数与目标分布输出更新后的参数 [q_mu, q_sigma2] Func_CVB_biGauss(mu_true, Sigma_true, q_mu, q_sigma2, 0.5);这里mu_true和Sigma_true是我们要逼近的目标对应后验分布的真实参数q_mu和q_sigma2是变分分布自身的均值与边缘方差最后传入的0.5是更新步长或者学习率一类的松弛系数。sigma2SIGMA.m这个辅助函数的存在感很强——它负责把一维的边缘方差还原成完整的协方差矩阵说明 CVB 的双高斯实现是显式地在维护“边缘分布参数”和“依赖结构参数”这两套量。2.3 混合高斯场景与蒙特卡洛评估流程双高斯场景跑通之后再来看混合高斯场景。它的流程比双高斯复杂因为引入了离散的聚类分配变量。主入口是A0_MAIN.m它内部会依次调用B0_SETTING.m设置模型参数用蒙特卡洛方式生成模拟数据然后让 CVB 和标准 VB 在同样的数据上分别做推断最后做对比评估。% 切换到混合高斯场景 cd(你的本地路径/Copula-Variational-Bayes-master/CVB for Gaussian mixture); % 运行混合高斯场景总入口 A0_MAIN这个场景里最值得关注的是c2_MonteCarlo.m。它做的事情是从 ground truth 分布里生成一批蒙特卡洛样本分别用 CVB 和 mean-field VB 做聚类推断再用聚类评估指标量化两个方法的差异。如果想复现论文里“CVB 在高相关成分下优于 mean-field VB”的结论重点看这个脚本里对数据生成的相关性参数是怎么设置的。2.4 怎么判断跑成功了判断运行是否成功不要只看有没有报错要看三个信号。第一命令行窗口有没有输出迭代过程中的 ELBO 或 KLD 数值且数值是否在递减或收敛到稳定值第二Figures目录下有没有生成对比图图里 CVB 的等高线是否比 VB 的等高线更贴近真实后验第三如果改了参数KLD 数值量级是否合理二维高斯之间的 KL 散度过大或为负都说明出了问题。3. 核心原理copula 如何把相关性塞回变分后验以及 Func_CVB_biGauss 的关键更新3.1 从 mean-field 到 copulaELBO 里多出的一项要理解 CVB 和普通 VB 的本质差异先看 mean-field 做了什么。标准变分贝叶斯假设后验分布可以因子分解为q(θ) ∏ qᵢ(θᵢ)也就是说所有参数在近似分布里是相互独立的。这个假设直接导致 ELBO 最大化问题退化成对每个参数独立做坐标上升计算变得简单但代价是丢失了参数之间的全部相关性。当真实后验里存在强相关时mean-field 逼近的结果会过度自信地低估不确定性而且均值估计本身也可能偏移。CVB 的切入点非常直接既然问题出在独立假设上那就用一个 copula 密度函数把各个边缘分布“绑”回联合分布q(θ) c(u₁, u₂, …, u_d) · ∏ qᵢ(θᵢ)其中 uᵢ Fᵢ(θᵢ) 是第 i 个边缘分布的累积概率变换c 是 copula 密度。这样改完之后ELBO 的展开式里会多出关键的一项% CVB 的 ELBO 分解结构示意 % log p(x) E_q[log p(x, theta)] - E_q[log q(theta)] % % 将 q(theta) c(u1,...,ud) * prod_i q_i(theta_i) 代入展开 % % ELBO E_q[log p(x, theta)] % - sum_i E_q[log q_i(theta_i)] % 边缘熵项 % - E_q[log c(u1,...,ud)] % copula 熵项 % % 与 mean-field 相比多出来的就是最后一项 % copula 熵项的梯度会把参数间的相关性反馈进每个边缘分布的更新这项的意义在于它把“分布之间的依赖结构”显式地变成了目标函数的一部分。当真实后验具有显著相关性时这一项的存在会引导变分分布往保留相关性的方向走这正是 CVB 相比 mean-field 的核心优势。而 information geometry 在这里的作用是为 copula 族的选择和参数更新提供几何视角——在概率分布流形上KL 散度和 Fisher 信息度量给了更新方向一个自然的定义。3.2 双高斯场景解析更新为什么能直接写出来双高斯场景是 CVB 最理想化的案例因为所有环节都存在解析表达式。真实后验是二维高斯变分分布在 copula 的帮助下也是二维高斯边缘是高斯高斯 copula 也是标准形式。于是 ELBO 里的每一项都能写成 μ 和 Σ 的代数式迭代更新就是普通的矩阵运算。下面这段代码是我按照工程里Func_CVB_biGauss.m思路整理的更新结构骨架方便理解每一行在做什么% 双高斯 CVB 单步更新示意结构具体数值以工程内实现为准 function [mu_new, sigma2_new] cvb_bigauss_update(mu0, Sigma0, mu_q, sigma2_q, rho) % mu0, Sigma0 : 目标分布或先验的均值与协方差 % mu_q, sigma2_q : 当前变分分布的均值与边缘方差 % rho : 变分分布内部的相关系数由 copula 参数决定 % 第一步由边缘方差和相关系数拼出变分协方差矩阵 Sigma_q [sigma2_q(1) rho*sqrt(sigma2_q(1)*sigma2_q(2)); rho*sqrt(sigma2_q(1)*sigma2_q(2)) sigma2_q(2)]; % 第二步精度矩阵相加对应信息形式的 VB 更新 Lambda_new inv(Sigma0) inv(Sigma_q); % 第三步均值按精度加权合并 mu_new Lambda_new \ (inv(Sigma0)*mu0(:) inv(Sigma_q)*mu_q(:)); % 第四步更新后的协方差矩阵转回边缘方差 Sigma_new inv(Lambda_new); sigma2_new diag(Sigma_new).; end这段代码里最关键的是rho这个参数。在 mean-field VB 里变分分布内部相关系数被强制为 0对应Sigma_q退化为对角矩阵而在 CVB 里rho是自由参数由 copula 决定并且会随着迭代更新。换句话说mean-field 是 CVB 在 rho0 处的一个特例。工程里sigma2SIGMA.m函数做的正是从边缘方差到完整协方差矩阵的组装与这里第一步做的事情相同。参数方面值得注意mu0和Sigma0既可以理解为真实后验参数也可以理解为贝叶斯框架下的先验参数在双高斯场景中代码直接把它当作目标分布来逼近所以更新方向是让变分参数向它靠拢。Lambda_new是精度矩阵信息形式的高斯更新让代码实现非常干净。3.3 混合高斯场景Func_CVB123 里的分层更新混合高斯场景比双高斯复杂一个量级因为隐变量分成了两类离散的聚类分配变量 z 和连续的组件参数。Func_CVB123.m里这个名字里的 123我理解对应的是整个更新流程中分组处理的三类核心参数——组件均值、组件精度、以及分配变量的变分后验。实际以代码内部注释为准。B0_SETTING.m里集中了所有超参数包括混合分量数 K、先验均值、先验协方差尺度、迭代次数等。蒙特卡洛数据由c2a_generate_monte_data.m生成c2b_clusterplot.m负责把聚类结果可视化。整个流程里有一处容易被误读的设计Func_emMU.m和Func_emL.m两个文件名里带 em但这里并不是在跑 EM 算法而是用它们计算似然相关的量或矩估计辅助 CVB 的更新。看到函数名带 EM 就以为是 EM 算法会绕很大弯路。混合场景中 CVB 相对 mean-field 的优势比双高斯更明显当混合成分之间的分离度不高、数据点的真实归属存在模糊性时mean-field 对分配变量的独立假设会过度自信地给出接近 0 或 1 的后验概率而 CVB 通过 copula 保留了不同分配之间的弱相关性后验概率会更平滑、更接近真实值。c2_MonteCarlo.m里的重复实验对比本质上就是在量化这个差异。4. 不收敛、负 KLD、乱码注释CVB 调试中五个高频翻车点与排查顺序4.1 现象迭代几步后 KLD 发散甚至出现 NaNA0_MAIN_biGauss.m跑起来之后前几步 ELBO 还在下降到某个迭代突然跳成NaN或者 KLD 值直接飙到 1e10 这种量级。等高线图上 CVB 的椭圆越画越离谱。原因基本出在初始化上。q_sigma2初始值给到 0 附近或者rho超出 [-1,1] 的合法区间协方差矩阵在第一步组装时就接近奇异或非正定。后续更新涉及矩阵求逆非正定矩阵直接让结果崩掉。MATLAB 的inv不会主动报错而是给出一堆巨大数值。解决方法是初始化之后立刻检查协方差的正定性别等到迭代崩了才回头查% 初始化后立刻做正定性检查 Sigma_q sigma2SIGMA(q_sigma2); [~, flag] chol(Sigma_q); if flag ~ 0 error(初始协方差矩阵不正定请检查 q_sigma2 初始化或 rho 取值); end % 经验上q_sigma2 初值取目标方差的一半量级即可 % rho 初值取 0 附近比较稳比如 0.2迭代后再让它自己漂移4.2 现象两个多元高斯之间的 KLD 算出来是负的Func_KLDMultiGauss.m输出的 KL 散度是负数。KL 散度理论上永远大于等于 0所以第一反应应该是代码调用有问题而不是理论有问题。原因通常是计算式里三项的符号或者传参顺序发生了混乱。多元高斯 KL 散度的标准形式包含三项协方差迹项、均值二次型项、行列式对数项。常见错误是把tr(Σ₂⁻¹Σ₁)写成了tr(Σ₁⁻¹Σ₂)或者行列式对数项里两个矩阵的顺序颠倒导致结果整体偏移。解决的办法是先手动构造一个已知结果校验一下% 一维情况手动验证N(0,1) 对 N(0,4) % 理论值 0.5 * (log(4) - 1 1/4) 0.4431 KL_ref Func_KLDMultiGauss([0], 1, [0], 4); disp(KL_ref); % 应该约等于 0.4431如果这个值对不上说明函数内部公式的项次顺序需要核对。改的时候注意保持行列式项和迹项的一致性不要只改一项。4.3 现象在 MATLAB 2023b 里打开 .m 文件中文注释全部乱码我在工程里加了一些中文注释重新打开之后全部变成乱码。这个在新版 MATLAB 里非常常见原因是新版编辑器默认按 UTF-8 解析文件而旧编辑器或 Windows 记事本保存的文件通常是 GBK/GB2312 编码。文件里的中文注释被错误解码后直接显示成乱码。解决方法是统一编码格式。最简单的是在 MATLAB 里重新打开文件另存为 UTF-8也可以先用命令行工具批量转换# 在项目根目录下用 iconv 批量把 .m 文件从 GBK 转成 UTF-8 find . -name *.m -exec iconv -f GBK -t UTF-8 {} -o {}.tmp \; -exec mv {}.tmp {} \;不过这个工程本身的注释是英文的一般不会触发乱码。如果你只是自己加中文注释建议直接在 MATLAB 编辑器里把默认编码设为 UTF-8路径在预设项里的“编辑器/Debugger”相关设置中调整。从那以后我只要拿到一个老工程的 .m 文件第一件事就是确认它的保存编码。4.4 现象等高线图空白或者只画出后半截椭圆C0_PLOT.m或者D0_PLOT.m输出的图里轮廓线要么完全空白要么只画出一半边缘被切断。原因是绘图网格的范围没有盖住分布的真实支撑区域。等高线的计算依赖meshgrid生成的网格如果网格只覆盖了坐标平面上很小一块区域而真实后验的协方差较大、分布半径超出网格边界画出来的就只是分布的一个局部切片。解决方法是先打印目标分布的范围再根据均值加减若干倍标准差设定网格边界% 先看目标分布的范围 fprintf(mu_true [%g, %g]\n, mu_true(1), mu_true(2)); fprintf(Sigma_true 对角线 [%g, %g]\n, Sigma_true(1,1), Sigma_true(2,2)); % 网格范围均值 ± 3.5 倍标准差覆盖 99% 以上的概率质量 x_range mu_true(1) 3.5*sqrt(Sigma_true(1,1)) * [-1 1]; y_range mu_true(2) 3.5*sqrt(Sigma_true(2,2)) * [-1 1]; [X, Y] meshgrid(linspace(x_range(1), x_range(2), 200), ... linspace(y_range(1), y_range(2), 200));网格点数 200 是经验值覆盖范围正确时画出来的等高线是平滑的闭合椭圆。如果出现锯齿把 200 提高到 400。4.5 现象蒙特卡洛重复实验的结果和 README 对不上混合高斯场景里c2_MonteCarlo.m观测到的聚类准确率和 README 描述的差异很大有时甚至平均后比标准 VB 还差。原因是随机数状态不一致。蒙特卡洛实验对随机数生成顺序极其敏感不同版本 MATLAB 的默认随机数生成器算法存在差异同一个randn调用序列在不同版本下产生不同的数据。ground truth 参数由c1_para_groundtruth.m定义但生成c2a_generate_monte_data.m的数据时没有固定种子每次运行数据都在变结果自然波动很大。解决办法是在运行脚本前固定随机数种子并且显式指定生成器% 固定随机数种子并显式指定 Mersenne Twister rng(2024, twister);固定种子之后整个蒙特卡洛实验变成可重复的。如果固定了种子结果仍然对不上再检查c2_MonteCarlo.m里评估指标的计算方式是否和 README 一致。5. 把 CVB 搬到自己的数据改维度、换边缘分布与三条验证手段5.1 从二维改到三维最小改动清单双高斯场景从二维扩展到三维需要改动的不是循环逻辑而是几处参数维度。mu_true从长度 2 变成 3Sigma_true变成 3x3q_sigma2从长度 2 变成 3。最隐蔽的坑是 copula 相关系数矩阵R二维时只有一个元素rho三维时有三个成对相关还要保证正定性。我的做法是先用下三角参数化构造相关矩阵绕开正定性问题% 构造三维高斯 copula 相关矩阵用 cholesky 因子保证正定 L [1 0 0; 0.3 0.95 0; 0.2 0.1 sqrt(1-0.2^2-0.1^2)]; R L * L.; % 归一化到对角为 1 R R ./ sqrt(diag(R) * diag(R).);这样构造出来的R天然满足正定条件不会出现第 4.1 节那种初始化就炸的问题。双高斯场景里rho是被当作自由参数迭代更新的但到了三维以上直接用 cvx 一类的工具重新做投影也行如果追求简单固定R只更新边缘参数CVB 仍然比 mean-field 多出一截表达能力。5.2 换边缘分布从高斯到 t 分布的注意点如果要把边缘分布从高斯换成 t 分布ELBO 里的熵项不再有解析解这是最大的变化。高斯场景里那些干净的矩阵更新式全部失效需要转用数值积分近似 ELBO 里的期望项。常见做法是用 Gauss-Hermite 求积来估计一维积分% 用 Gauss-Hermite 近似计算 t 分布边缘的熵示意 % 节点和权重可以预先算好比如 32 点 [xi, wi] gauss_hermite(32); % 熵近似E[log q(theta)] ≈ sum_i wi .* log(tpdf(mu sqrt(2)*nu*xi, dof)) % 其中 nu 是 t 分布的尺度参数 entropy_q sum(wi .* log(tpdf(mu_q sqrt(2)*nu*xi, dof)));一维用 Gauss-Hermite 精度还不错但注意它假设了高斯权重t 分布尾部较重时节点覆盖范围可能不够需要把节点乘上一个大于 1 的膨胀因子再试。换分布之后KLD 的计算也必须换成对应的数值版本不能再复用Func_KLDMultiGauss.m。5.3 三条验证手段如何判断近似质量换数据换分布之后怎么知道 CVB 的结果值不值得信任我一般强制自己走三条验证缺一不可。第一条是 KLD 数值对照。算一下真实分布与 CVB 近似的 KLD再算一下真实分布与 mean-field 近似的 KLD看 CVB 是否确实更小。如果两种方法差距在 5% 以内说明这个场景下相关性并不强CVB 的收益有限。第二条是等高线视觉对照。把真实后验和两种近似的等高线画在同一张图里重点看椭圆的长轴方向是否对齐。KLD 数值可能有误差但视觉不会骗人长轴方向对不上就是依赖结构没学对。第三条是边缘分布验证。把变分分布的边缘 CDF 和真实后验的边缘 CDF 做 QQ plot如果点在直线上说明边缘近似得好。这条专门防止 CVB 出现“联合分布对但边缘偏了”的畸形结果。Copula 建模的初衷就是边缘和依赖结构分开处理边缘错了 copula 再准也白搭。CVB 是一个比 mean-field 更接近真实后验的框架但工程落地时那些矩阵正定、编码格式、随机种子的小问题能靠一个下午的调试轻松消磨掉你的耐心。我从第一次被负 KLD 坑过之后现在拿到任何新数据都会强制先跑一遍 KLD 对照、等高线对比和边缘 QQ plot 这三条验证再决定要不要信这个结果。希望帮到你。本文还有配套的精品资源点击获取