wehub-resource-sync
78 lines
1.7 KiB
Matlab
Executable File
78 lines
1.7 KiB
Matlab
Executable File
function [model, energy] = linRegVb(X, t, prior)
|
|
% Variational Bayesian inference for linear regression.
|
|
% Input:
|
|
% X: d x n data
|
|
% t: 1 x n response
|
|
% prior: prior parameter
|
|
% Output:
|
|
% model: trained model structure
|
|
% energy: variational lower bound
|
|
% Written by Mo Chen (sth4nth@gmail.com).
|
|
if nargin < 3
|
|
a0 = 1e-4;
|
|
b0 = 1e-4;
|
|
c0 = 1e-4;
|
|
d0 = 1e-4;
|
|
else
|
|
a0 = prior.a;
|
|
b0 = prior.b;
|
|
c0 = prior.c;
|
|
d0 = prior.d;
|
|
end
|
|
[m,n] = size(X);
|
|
I = eye(m);
|
|
xbar = mean(X,2);
|
|
tbar = mean(t,2);
|
|
|
|
X = bsxfun(@minus,X,xbar);
|
|
t = bsxfun(@minus,t,tbar);
|
|
|
|
XX = X*X';
|
|
Xt = X*t';
|
|
|
|
maxiter = 100;
|
|
energy = -inf(1,maxiter+1);
|
|
tol = 1e-8;
|
|
|
|
a = a0+m/2; % 10.94
|
|
c = c0+n/2;
|
|
Ealpha = 1e-4;
|
|
Ebeta = 1e-4;
|
|
for iter = 2:maxiter
|
|
% q(w)
|
|
invS = diag(Ealpha)+Ebeta*XX; % 10.101
|
|
U = chol(invS);
|
|
Ew = Ebeta*(U\(U'\Xt)); % 10.100
|
|
KLw = -sum(log(diag(U)));
|
|
% q(alpha)
|
|
w2 = dot(Ew,Ew);
|
|
invU = U'\I;
|
|
trS = dot(invU(:),invU(:));
|
|
b = b0+0.5*(w2+trS); % 10.95
|
|
Ealpha = a/b; % 10.102
|
|
KLalpha = -a*log(b);
|
|
% q(beta)
|
|
e2 = sum((t-Ew'*X).^2);
|
|
invUX = U'\X;
|
|
trXSX = dot(invUX(:),invUX(:));
|
|
d = d0+0.5*(e2+trXSX);
|
|
Ebeta = c/d;
|
|
KLbeta = -c*log(d);
|
|
% lower bound
|
|
energy(iter) = KLalpha+KLbeta+KLw;
|
|
if energy(iter)-energy(iter-1) < tol*abs(energy(iter-1)); break; end
|
|
end
|
|
const = gammaln(a)-gammaln(a0)+gammaln(c)-gammaln(c0)+a0*log(b0)+c0*log(d0)+0.5*(m-n*log(2*pi));
|
|
energy = energy(2:iter)+const;
|
|
w0 = tbar-dot(Ew,xbar);
|
|
|
|
model.w0 = w0;
|
|
model.w = Ew;
|
|
model.alpha = Ealpha;
|
|
model.beta = Ebeta;
|
|
model.a = a;
|
|
model.b = b;
|
|
model.c = c;
|
|
model.d = d;
|
|
model.xbar = xbar;
|