My favorites | Sign in
Project Home
READ-ONLY: This project has been archived. For more information see this post.
Search
for
  Advanced search   Search tips   Subscriptions
Issue 4: SavageDickey with several sources
1 person starred this issue and may be notified of changes. Back to list
Status:  New
Owner:  ----


 
Reported by lopez.al...@gmail.com, Nov 26, 2013
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
Project Member #1 jean.dau...@gmail.com
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.


Dec 3, 2013
#2 lopez.al...@gmail.com
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
Project Member #3 lionel.r...@gmail.com
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
#4 lopez.al...@gmail.com
Thank you very much Lionel. 

It seems that works !

Best, 

Alizée

Powered by Google Project Hosting