forked from translationalneuromodeling/tapas
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtapas_h2gf_prepare_ptheta.m
142 lines (118 loc) · 4.27 KB
/
tapas_h2gf_prepare_ptheta.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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
function [ptheta] = tapas_h2gf_prepare_ptheta(ptheta, theta, pars)
%%
%
% copyright (C) 2016
%
% Check input
[nr, nc] = size(ptheta.c_prc.priormus);
assert(nc == 1 || nr == 1, 'tapas:h2gf:input', ...
'Priors should be vectors');
if nc > 1
ptheta.c_prc.priormus = ptheta.c_prc.priormus';
end
[nr, nc] = size(ptheta.c_prc.priorsas);
assert(nc == 1 || nr == 1, 'tapas:h2gf:input', ...
'Priors should be vectors');
if nc > 1
ptheta.c_prc.priorsas = ptheta.c_prc.priorsas';
end
[nr, nc] = size(ptheta.c_obs.priormus);
assert(nc == 1 || nr == 1, 'tapas:h2gf:input', ...
'Priors should be vectors');
if nc > 1
ptheta.c_obs.priormus = ptheta.c_obs.priormus';
end
[nr, nc] = size(ptheta.c_obs.priorsas);
assert(nc == 1 || nr == 1, 'tapas:h2gf:input', ...
'Priors should be vectors');
if nc > 1
ptheta.c_obs.priorsas = ptheta.c_obs.priorsas';
end
% Number of perceptual parameters
i = 1;
n = numel(ptheta.c_prc.priormus);
% Indices of perceptual parameters
ptheta.theta_prc = i : n + i -1;
% Number of perceptual parameters
ptheta.theta_prc_nd = n;
% Number of observation parameters
i = i + n;
n = numel(ptheta.c_obs.priormus);
% Indices of observation parameters
ptheta.theta_obs = i : i + n - 1;
% Number of observation parameters
ptheta.theta_obs_nd = n;
% Prior means
ptheta.mu = [ptheta.c_prc.priormus; ptheta.c_obs.priormus];
% Prior precisions
ptheta.pe = 1./[ptheta.c_prc.priorsas; ptheta.c_obs.priorsas];
% Initial sampling point
ptheta.p0 = ptheta.mu;
% Indices of sampled (ie, non-fixed) perceptual parameters
prc_idx = ptheta.c_prc.priorsas;
prc_idx(isnan(prc_idx)) = 0;
prc_idx = find(prc_idx);
% Indices of sampled (ie, non-fixed) observation parameters
obs_idx = ptheta.c_obs.priorsas;
obs_idx(isnan(obs_idx)) = 0;
obs_idx = find(obs_idx);
% Indices of all sampled (ie, non-fixed) parameters
valid = [prc_idx; numel(ptheta.c_prc.priorsas) + obs_idx];
% Projection matrix from all-parameter to sampled-parameter space
ptheta.jm = zeros(size(ptheta.mu, 1), size(valid, 1));
for i = 1:numel(valid)
ptheta.jm(valid(i), i) = 1;
end
ptheta.jm = sparse(ptheta.jm);
% Remove from p0 the prior in all active rows. For the non active rows,
% This line does not have an effect.
ptheta.p0 = ptheta.p0 - (ptheta.jm * ptheta.jm') * ptheta.p0;
% Stored a lower dimensional representation with only the valid or active
% parameters
ptheta.mu = ptheta.jm' * ptheta.mu;
ptheta.pe = ptheta.jm' * ptheta.pe;
% Check for empirical prior
if ~isfield(ptheta, 'empirical_priors')
ptheta.empirical_priors = struct();
end
% Check for weighting of the prior
if ~isfield(ptheta.empirical_priors, 'eta')
% Defaults to 1
ptheta.empirical_priors.eta = ptheta.jm' * ones(size(ptheta.jm, 1), 1);
elseif isscalar(ptheta.empirical_priors.eta)
% If eta is a scalar, use the same eta for all parameters.
ptheta.empirical_priors.eta = ptheta.jm' * ...
ptheta.empirical_priors.eta * ...
ones(size(ptheta.jm, 1), 1);
else
% Otherwise it should be a vector of dimensions perceptual + observational
ptheta.empirical_priors.eta = ptheta.jm' * ptheta.empirical_priors.eta;
end
% Hyperprior of the population mean
if ~isfield(ptheta.empirical_priors, 'alpha')
ptheta.empirical_priors.alpha = (ptheta.empirical_priors.eta + 1)./2;
elseif isscalar(ptheta.empirical_priors.alpha)
% If eta is a scalar, use the same eta for all parameters.
ptheta.empirical_priors.eta = ptheta.jm' * ...
ptheta.empirical_priors.alpha * ...
ones(size(ptheta.jm, 1), 1);
else
% Otherwise it should be a vector of dimensions perceptual + observational
ptheta.empirical_priors.alpha = ptheta.jm' * ptheta.empirical_priors.alpha;
end
if ~isfield(ptheta.empirical_priors, 'beta')
% Get the mean of the prior provided by the user.
pe = ptheta.pe;
% Defaults to 1
ptheta.empirical_priors.beta = ptheta.empirical_priors.eta ./ (2.0 * pe);
elseif isscalar(ptheta.empirical_priors.beta)
% If beta is a scalar, use the same beta for all parameters.
ptheta.empirical_priors.beta = ptheta.jm' * ...
ptheta.empirical_priors.beta * ...
ones(size(ptheta.jm, 1), 1);
else
% Otherwise it should be a vector of dimensions perceptual + observational
ptheta.empirical_priors.beta = ptheta.jm' * ptheta.empirical_priors.beta;
end
end