SavageDickey with several sources

[lionel-rigoux]

From @GoogleCodeExporter on October 7, 2015 9:29

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

Original issue reported on code.google.com by lopez.al...@gmail.com on 26 Nov 2013 at 1:42

Copied from original issue: lionel-rigoux/mbb-vb-toolbox#4

[lionel-rigoux]

From @GoogleCodeExporter on October 7, 2015 9:29

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.

Original comment by jean.dau...@gmail.com on 3 Dec 2013 at 9:48

[lionel-rigoux]

From @GoogleCodeExporter on October 7, 2015 9:29

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.

Original comment by lopez.al...@gmail.com on 3 Dec 2013 at 11:08

[lionel-rigoux]

From @GoogleCodeExporter on October 7, 2015 9:29

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

Original comment by lionel.r...@gmail.com on 4 Dec 2013 at 2:53

[lionel-rigoux]

From @GoogleCodeExporter on October 7, 2015 9:29

Thank you very much Lionel. 

It seems that works !

Best, 

Alizée

Original comment by lopez.al...@gmail.com on 4 Dec 2013 at 6:32