-
Notifications
You must be signed in to change notification settings - Fork 0
/
spm_BMS_F.m
44 lines (37 loc) · 1.2 KB
/
spm_BMS_F.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
function [F_samp,F_bound] = spm_BMS_F(alpha,lme,alpha0)
% Compute two lower bounds on model evidence p(y|r) for group BMS
% FORMAT [F_samp,F_bound] = spm_BMS_F(alpha,lme,alpha0)
%
% alpha - parameters of p(r|y)
% lme - array of log model evidences
% rows: subjects
% columns: models (1..Nk)
% alpha0 - priors of p(r)
%
% F_samp - sampling estimate of <ln p(y_n|r)>
% F_bound - lower bound on lower bound of <ln p(y_n|r)>
%
% Reference:
% Stephan KE, Penny WD, Daunizeau J, Moran RJ, Friston KJ
% Bayesian Model Selection for Group Studies. Neuroimage 2009 46(4):1004-17
%__________________________________________________________________________
% Will Penny
% Copyright (C) 2008-2022 Wellcome Centre for Human Neuroimaging
alpha0 = sort(alpha0);
if alpha0(1) ~= alpha0(end)
error('alpha0 should have identical values.')
end
alpha0 = alpha0(1);
a_sum = sum(alpha);
psi_sum = psi(a_sum);
psi_diff = psi(alpha) - psi_sum;
gm = gammaln(alpha);
[s_samp,s_bound] = spm_BMS_F_smpl(alpha,lme,alpha0);
K = length(alpha);
F = 0;
for k = 1:K
F = F - (alpha(k) - alpha0)*psi_diff(k) + gm(k);
end
F = F - gammaln(a_sum);
F_bound = F + s_bound;
F_samp = F + s_samp;