-
Notifications
You must be signed in to change notification settings - Fork 10
/
tlsa_decode_gaussian.m
71 lines (62 loc) · 2.24 KB
/
tlsa_decode_gaussian.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
function [mu iS] = tlsa_decode_gaussian(data,testdata,results,mu0,iS0)
% Decode continuous covariates from brain images under Gaussian assumptions.
% If no prior parameters are provided, an empirical prior is constructed using
% the sample mean and precision from the training data.
%
% USAGE: [mu iS] = tlsa_decode_gaussian(data,testdata,results,[mu0],[iS0])
%
% INPUTS:
% data - [1 x S] training data structure
% testdata - [1 x S] test data structure
% results - results structure
% mu0 (optional) - prior mean
% iS0 (optional) - prior precision matrix
%
% OUTPUTS:
% mu - [1 x S] cell array, where each cell contains the [N x C] posterior
% mean covariate vector for each datapoint
% iS - [1 x S] cell array, where each cell contains the [C x C] posterior precision matrix
%
% Sam Gershman, Oct 2012
for s = 1:length(data)
% if no mean & precision given, estimate empirical (MLE) prior
if nargin < 4
[mu0 iS0] = tlsa_empirical_prior(data(s));
end
q = results.q(s);
tau = q.rho*q.nu;
F = tlsa_map(results.opts.mapfun,q.omega,data(s).R);
A = q.W*F;
iS{s} = iS0 + tau*(A*A');
mu{s} = bsxfun(@plus,mu0*iS0,tau*testdata(s).Y*A')/iS{s};
end
end
function [mu iS] = tlsa_empirical_prior(data,diagonalize)
% Constructs an empirical prior for the covariates based on the training data.
%
% This prior is equivalent to the maximum likelihood estimates of
% the mean and covariance matrix of a Gaussian.
%
% USAGE: [mu iS] = tlsa_empirical_prior(data,[diagonalize])
%
% INPUTS:
% data - data structure
% diagonalize (optional) - if 1, uses diagonal approximation of the
% precision matrix (default: 1)
%
% OUTPUTS:
% mu - [1 x C] mean vector
% iS - [C x C] precision matrix
%
% Sam Gershman, Sep 2011
if nargin < 2 || isempty(diagonalize); diagonalize = 1; end
X = data.X;
mu = mean(X);
r = bsxfun(@minus,X,mu);
Sigma = (r'*r)/size(X,1);
if diagonalize
iS = diag(1./(eps+diag(Sigma)));
else
iS = inv(Sigma);
end
end