| Issue 4: | SavageDickey with several sources | |
| 1 person starred this issue and may be notified of changes. | Back to list |
What steps will reproduce the problem?
1. VBA_NLStateSpaceModel with options.sources is a 1x3 struct array (two continuous data and one binary data)
2. VB_SavageDickey on the result of the first inversion
What is the expected output? What do you see instead?
Output expected : free energy and approximate log evidence of my reduced model.
Error seen :
Error using ^
Inputs must be a scalar and a square matrix.
To compute elementwise POWER, use POWER (.^) instead.
Error in VB_SavageDickey (line 68)
Sf = po1.a_sigma./(po1.b_sigma^2);
Error in test_param_all_tasks (line 27)
[F2(s,h),po2(s,h)]=VB_SavageDickey(po1,pr1,F1,dim1,pr2);
What version of the product are you using? On what operating system?
Version used : 1074 (last commit revision)
Operating system : Windows
Matlab version : 2013a
Please provide any additional information below.
If I modify the file VB_SavageDickey at line 68, 70 and 72 by remplacing "^" by ".^", I have a new error :
Index exceeds matrix dimensions.
Error in spm_log_evidence (line 92)
qP = spm_inv(qC(i,i),TOL);
Error in VB_SavageDickey (line 73)
[dF,mr,Sr] = spm_log_evidence(mf,Sf,mf0,Sf0,mr0,Sr0);
Error in test_param_all_tasks (line 27)
[F2(s,h),po2(s,h)]=VB_SavageDickey(po1,pr1,F1,dim1,pr2);
Is there an other way to use VB_SavageDickey with several sources ?
Thanks in advance,
Alizée
Dec 3, 2013
Hey Jean,
Thank you very much for your answer.
I replaced the section you indicated and I have a new error, here it is :
%%%%%%%
Index exceeds matrix dimensions.
Error in spm_log_evidence (line 92)
qP = spm_inv(qC(i,i),TOL);
Error in VB_SavageDickey (line 82)
[dF,mr,Sr] = spm_log_evidence(mf,Sf,mf0,Sf0,mr0,Sr0);
Error in test_param_all_tasks (line 27)
[F2(s,h),po2(s,h)]=VB_SavageDickey(po1,pr1,F1,dim1,pr2);
%%%%%%
It seems it is still in the same section, any idea ?
Cheers,
Alizée.
Dec 4, 2013
It seems that the spm function cannot deal with vector arguments. One solution is to loop over sources: replace the code
%
[dF,mr,Sr] = spm_log_evidence(mf,Sf,mf0,Sf0,mr0,Sr0);
po2.b_sigma = mr./Sr;
po2.a_sigma = po2.b_sigma.*mr;
F2 = F2 +dF;
%
by the following loop:
%
for iSource = 1:numel(po1.a_alpha)
[dF,mr,Sr] = spm_log_evidence(mf(iSource), ...
Sf(iSource),mf0(iSource),Sf0(iSource),mr0(iSource),Sr0(iSource));
po2.b_alpha(iSource) = mr/Sr;
po2.a_alpha(iSource) = po2.b_alpha(iSource)*mr;
F2 = F2 +dF;
end
%
This should do the trick! Tell us if it works so we can commit the changes.
Best,
Lionel
Dec 4, 2013
Thank you very much Lionel. It seems that works ! Best, Alizée |
Hey Alizee, This issue is due to the fact that some advanced functionalities (like VB_SavageDickey) were not modified to account for multi-source inversions. More precisely, multi-source inversion can be used with priors on noise precision that are the same for all sources, but results in source-specific posteriors. This means that the size of the priors and the posterior may not be the same! The following patch should work: --> replace the section starting with 'if isfield(po1,...' (line 66 in VB_SavageDickey.m) with the following piece of code: if isfield(po1,'a_sigma') && ~isempty(po1.a_sigma) ns = length(po1.b_sigma); mf = po1.a_sigma./po1.b_sigma; Sf = diag(po1.a_sigma./(po1.b_sigma.^2)); mf0 = pr1.a_sigma./pr1.b_sigma; Sf0 = pr1.a_sigma./(pr1.b_sigma.^2); if ~isequal(length(mf0),ns) mf0 = mf0.*ones(ns,1); Sf0 = Sf0*eye(ns); end mr0 = pr2.a_sigma./pr2.b_sigma; Sr0 = pr2.a_sigma./(pr2.b_sigma.^2); if ~isequal(length(mr0),ns) mr0 = mr0.*ones(ns,1); Sr0 = Sr0*eye(ns); end [dF,mr,Sr] = spm_log_evidence(mf,Sf,mf0,Sf0,mr0,Sr0); po2.b_sigma = mr./Sr; po2.a_sigma = po2.b_sigma.*mr; F2 = F2 +dF; end Tell me whetehr this solves your problem. If it does, I will commit the change... Cheers, Jean.