wehub-resource-sync
60 lines
1.3 KiB
Matlab
Executable File
60 lines
1.3 KiB
Matlab
Executable File
function [model, llh] = linRegEm(X, t, alpha, beta)
|
|
% Fit empirical Bayesian linear regression model with EM (p.448 chapter 9.3.4)
|
|
% Input:
|
|
% X: d x n data
|
|
% t: 1 x n response
|
|
% alpha: prior parameter
|
|
% beta: prior parameter
|
|
% Output:
|
|
% model: trained model structure
|
|
% llh: loglikelihood
|
|
% Written by Mo Chen (sth4nth@gmail.com).
|
|
if nargin < 3
|
|
alpha = 0.02;
|
|
beta = 0.5;
|
|
end
|
|
[d,n] = size(X);
|
|
I = eye(d);
|
|
xbar = mean(X,2);
|
|
tbar = mean(t,2);
|
|
|
|
X = bsxfun(@minus,X,xbar);
|
|
t = bsxfun(@minus,t,tbar);
|
|
|
|
XX = X*X';
|
|
Xt = X*t';
|
|
|
|
tol = 1e-4;
|
|
maxiter = 100;
|
|
llh = -inf(1,maxiter+1);
|
|
for iter = 2:maxiter
|
|
A = beta*XX+alpha*eye(d);
|
|
U = chol(A);
|
|
|
|
m = beta*(U\(U'\Xt));
|
|
m2 = dot(m,m);
|
|
e2 = sum((t-m'*X).^2);
|
|
|
|
logdetA = 2*sum(log(diag(U)));
|
|
llh(iter) = 0.5*(d*log(alpha)+n*log(beta)-alpha*m2-beta*e2-logdetA-n*log(2*pi)); % 3.86
|
|
if abs(llh(iter)-llh(iter-1)) < tol*abs(llh(iter-1)); break; end
|
|
|
|
invU = U'\I;
|
|
trS = dot(invU(:),invU(:)); % A=inv(S)
|
|
alpha = d/(m2+trS); % 9.63
|
|
|
|
invUX = U'\X;
|
|
trXSX = dot(invUX(:),invUX(:));
|
|
beta = n/(e2+trXSX); % 9.68 is wrong
|
|
end
|
|
w0 = tbar-dot(m,xbar);
|
|
|
|
llh = llh(2:iter);
|
|
model.w0 = w0;
|
|
model.w = m;
|
|
%% optional for bayesian probabilistic inference purpose
|
|
model.alpha = alpha;
|
|
model.beta = beta;
|
|
model.xbar = xbar;
|
|
model.U = U;
|