变分高斯 Copula 推理（Matlab代码实现）

    目录

💥1 概述

📚2 运行结果

🎉3 参考文献

👨‍💻4 Matlab代码

💥1 概述

​变分推理是计算后验分布的一种常用方法,因其收敛速度快、理论基础扎实,在机器学习中发挥着核心作用。传统变分推理假设变分分布为平均场形式(Mean-field),即隐变量之间相互独立,这种理想化假设虽然可以简化优化过程、提高计算可行性,却忽略了隐变量之间的后验相关性(后验依赖关系),对具有变量相关性的模型来说并不适用。Copula变分推理(Copula Variational Inference,CVI)是解决上述问题的主流方法之一,通过使用Copula函数来捕获模型中隐变量之间的依赖关系,虽然可以获得全面且完整的信息,但存在以下问题:(1)由于Copula函数未考虑到依赖关系的稀疏性,从而捕获到一些非必要依赖关系,降低了变分推理近似的准确性;(2)Copula函数采用满秩的方法使得计算复杂度为隐变量数量的平方,增加了计算复杂性。本文针对以上两个问题对CVI方法进行研究和改进,提出了稀疏Copula变分推理方法,又因该方法对稀疏后的相关性信息获取不够全面,进而提出了稀疏Copula层次化变分推理方法。

本文使用高斯 copulas（结合固定/自由形式的边距）作为自动推理引擎，用于通用分层贝叶斯模型中的变分逼近（仅有的两个特定于模型的术语是对数似然和先验项及其导数）。本文还评估了在单变量边缘中再现的特殊性以及跨潜在变量广泛捕获的后验依赖性。

📚2 运行结果

主函数部分代码：

%%  VGC Numeric
opt.adaptivePhi = 0;  opt.normalize = 1;
fix.c = 1;  %
 unnormalizing constantopt.k = 10; % Degree/Maximum Degree of Bernstein Polynomials
opt.MaxIter = 50; % Number of SGD stages
opt.NumberZ = 1; % Average Gradients
opt.InnerIter = 500; % Number of iteration
opt.N_mc = 1; % Number of median average in sELBO stepsize search
opt.nsample = 1e5;
​
BPtype = 'BP';  % Bernstein Polynomials
% BPtype = 'exBP';  % Extended Bernstein Polynomials
PsiType = 'Normal';   opt.PsiPar(1) = 0;  opt.PsiPar(2) = 1;  % variance
​
opt.PoM.PsiType1 = 'Normal';
opt.PoM.PsiPar1(1) = 0;
opt.PoM.PsiPar1(2) = 1;  % variance
opt.PoM.PsiType2 = 'Exp';
opt.PoM.PsiPar2 = 1;
​
% PsiType = 'Exp';    opt.PsiPar = 1;
PhiType = 'Normal';  opt.PhiPar(1) = 0; opt.PhiPar(2) =1;  % variance
opt.Wthreshold = 1e4;
​
% Learning rate
opt.LearnRate.Mu = 1e-3; opt.LearnRate.C = 1e-4;
opt.LearnRate.W = 0.25.*1e-3; opt.LearnRate.dec = 0.95; % decreasing base learning rate
​
switch BPtype
    case 'BP'
        opt.D = opt.k;
    case 'exBP'
        opt.D = opt.k*(opt.k+1)/2; % # of basis functions
end
​
% Diagonal constraint on Upsilon
opt.diagUpsilon = 0;
​
%% Initialization
ini.Mu = opt.PhiPar(1).*ones(fix.P,1);
ini.C = sqrt(opt.PhiPar(2)).*0.1.*eye(fix.P);
%
 ini.w = randBPw(fix.P, opt.D, 
1
, 
1
);ini.w = 1./opt.D.*ones(fix.P, opt.D);
% ini.w = randBPw(fix.P, opt.D, 1, 1);
​
% Median Outlier Removal
opt.OutlierTol = 100; % Threshold for online outlier detection
opt.WinSize = 20; % Size of the window
ini.WinSet = ini_WinSet(trueModel, PhiType, PsiType, BPtype, fix, opt, ini);
​
%% VGC MuCW-Numeric
opt.updateW = 1; opt.updateMuC = 1;
[ELBO3, par3] = vgcbp_MuCW(trueModel, inferModel2, PhiType, PsiType, BPtype, VGmethod, fix, ini, opt);
​
tic;
figplot.Y_svc3 = sampleGC(PhiType, PsiType, BPtype, opt, par3);
figplot.Y_svc3(sum(~isfinite(figplot.Y_svc3),2)~=0,:) = [];
toc;
​
PoissonLogVGC = figplot.Y_svc3';
save PoissonLogVGC.mat PoissonLogVGC
%
% ------------------------------load PoissonLogVGC.mat
load PoissonLogMCMC.mat
​
beta0 = beta0(:)'; beta1 = beta1(:)';
beta2 = beta2(:)'; tau = tau(:)';
​
X_mcmc = [beta0; beta1; beta2; tau];
X_vgc =PoissonLogVGC;
ncoutour = 6; figplot.ls = 1.5; figplot.nbins = 50;
​
figure
for i = 1:4
    for j = 1:4
        k = (i-1)*4+j;
        subplot(4,4,k)
        if i==j
            [f_x1,x_x1] = hist(X_mcmc(i,:), figplot.nbins);
            [f_x2,x_x2] = hist(X_vgc(i,:), figplot.nbins);
            plot(x_x1, f_x1./trapz(x_x1,f_x1), 'r--', 'linewidth', 2)
            hold on
            plot(x_x2, f_x2./trapz(x_x2,f_x2), 'b-', 'linewidth', 2)
            hold off
            legend('JAGS','VGC-BP')
        else
            tmp = hist1D2D([X_mcmc(i,:); X_mcmc(j,:)]', figplot.nbins);
            tmp2 = hist1D2D([X_vgc(i,:); X_vgc(j,:)]'
, figplot.nbins);            contour(tmp.Seqx,tmp.Seqy, tmp.count./sum(sum(tmp.count)), ncoutour, 'r--' ,   'linewidth', 2*figplot.ls);
            hold on
            contour(tmp2.Seqx,tmp2.Seqy, tmp2.count./sum(sum(tmp2.count)), ncoutour, 'b' ,   'linewidth', 2*figplot.ls);
            hold off
            legend('JAGS', 'VGC-BP')
        end
    end
end

🎉3 参考文献

[1]曹竞月.Copula变分推理研究[D].吉林大学,2022.

部分理论引用网络文献，若有侵权联系博主删除。