mbb-team / vba-toolbox Goto Github PK

Dear developers,

My data is sampled irregularly ( frequently in the first several days, then the frequency dropped to once a week until once a month), is it possible for the VBA toolbox to deal with data like that?

If it is possible, could you please show how to set up the settings, please?

Thanks!

Yan

There are some redundancies in the core and display code to process differentially single or multi-sources data. Part of the switches between these two trails relies on the number of sources, some on an "extended" flag in the options structure. The latter was introduced to test the equivalence of the processing of single source models accross the two versions of codes. Therefore, this flag is not reliable to indicate the number of sources and only add unecessary redundancy in the instructions workflow.
I think the redundancy should be reduced to its minimum and the "extended" flag should be removed. In particular, this concerns the following functions:

• VBA_initDisplay.m / VBA_initDisplay_extended.m
• VBA_updateDisplay.m / VBA_updateDisplay_extended.m

Part of the redundancy is however required has the treatment of stochastic and multisources models cannot be yet fused (cf. VBA_Iphi.m)

Some spm functions rely on the nanxxx functions specific to spm. Those functions are also found in the statistics toolbox with a more general signature.
For the sake of independance, those functions should be included in the VBA toolbox.

MFX cannot be used to invert multisession data. As the multisession option is duplicating the priors, the posterior and prior do not have the same dimensions. This crashes MFX as it derives iteratively the prior from the posterior.
A potential solution would be to allow to specify priors for multisession data using the extended scheme (priors on all parameters, including session-wise duplicates).

In the documentation of the options for the priors it is said that iQy is the precision matrix of the observation error. As that is the inverse of the variance, then if the data scales with a factor of f, then iQy should scale by a factor of f^-2. All my experiments with the toolbox show however, that I need to scale iQy by f^-1in order to get equivalent results for scaled data, so for me it seems that iQy is rather a matrix square root of the precision matrix. Is that correct? If so, it would be nice to correct the documentation accordingly.

I try to perform multisession Volterra inversions. The way I currently do it in a non multisession way is to have a multidimensional y and put the inputs in inG. It will be quite qtraitforward for me to extend this inversion in a multisession way if it were possible to have one inG per session, basically to call inG{i} fot session i.
This would require a new field in the multisession structure to allow it if specified (default would be to use the same inG for all sessions).

For multiple reasons:

• the new display uses isYout, which is not passed correctly (ie. truncated) by onlinewrapper
• I broke the online display for hyperparameters (my bad)
• some suffstats are not passed correctly. In particular, as only one point/state is inverted, suffstats does not return any dx infos, only dx0.

I am creating a new branch fix-online to clean up the mess I induced by the change in display.
What is the expected behaviour? All params/HP are displayed as function of time?

In some cases the amount of clusters identified by RFT_main.m (nc line 165) and the number of peaks identified (is line 221) does not match. (If no local maxima could be found in particular because there is no peak in the cluster). Then RFT_main.m and hence also RFT_GLM_contrast.m crashes. Could this case be handled somehow, just by reporting there was no peak eventhough the cluster as a full is significant for example?

Hi VBA developers,
I recently had a chance to go through 'demo_VanDerPol' and 'f_vanDerPol' specifying the evolution function.
There are two things I want to discuss which may help improve the demo:

1. Practical assumption for the precision rates of hidden states

• Van Der Pol has two distinct hidden states (say x1 and x2), and the innovation for x1 is simply the replication of the second time series one time step earlier. The current demo assumes the precision rate (alpha) to be equal to both time series. I think that having the first time series being deterministic isn't a bad idea, by adding modified iQx{t}.

2. dF_dX in 'f_vanDerPol'

• This is more like a question, but I found that the Jacobian matrix for f(x) doesn't have a part for x in line 25-28. Is there reason that the partial derivatives of the evolution function only go on the innovation part? I was confused because in other scripts (e.g. line 18-19 in 'f_lin2D.m') I saw the entire f(x) partially derived with regard to x and had an identity matrix for the previous level. Can you clarify this?

I tried using the modified iQx{t} and dF_dX and got the better fit from the VBA inversion, hope this helps the consideration.

Thanks,
Jungmin

Multinomial data is not handled perfectly in the toolbox.

• axes of the predictive density are not correctly labelled
• demo and documentation for multinomial data sources are inexistant
• accuracy measures are missing or not documented

It would be good to add a line to check for NaN values in the VBA_orth function
Something like

if nonzeros( isnan(A) )
error('your matrix should not include any NaN value');
end

at the start of the script because otherwise any parasite NaN value will create a column full of NaN in the output possibly messing with other scripts afterwards.

I am using a non-informative prior setting for the gamma distributed precision of the dynamic noises in the second order differential model (VBA_StateSpaceModel). (priors.a_alpha =0.001; priors.b_alpha =0.001; priors.a_sigma =0.001; priors.b_sigma =0.001 ;) The mean value of the precision from the inferred results is in the order of 0.1 to 1, which is a quite big noise for my system.

Did I use a bad prior for the initial setting of the noise? How reliable is the inference of the hyperparameters for this toolbox?

Thank you in advance! I am looking forward to your reply.

Dear Lionel,

I installed VBA by running VBA_setup() on my windows pc in MATLAB R2022a but when running the Qlearning demo I get the following error

>> demo_Qlearning
No inputs provided, generating simulated behavior...

Simulating SDE... OK (took 0.129 seconds).
Deriving prior's sufficient statistics ... OK.
Main VB inversion...
VB iteration #1         F=-8.919e+01         ... dF=4.909e+00
Unrecognized function or variable 'nansum'.

Error in VBA_r2>sumall (line 51)
    s = nansum (nansum (z));

Error in VBA_r2 (line 43)
SS2_tot = sumall ((data(:) - VBA_nanmean(data(:))) .^2);

Error in VBA_getFit (line 88)
    fit.R2(si) = VBA_r2 (suffStat.gx(idx,:), out.y(idx,:), out.options.isYout(idx,:));

Error in VBA_wrapup (line 29)
out.fit = VBA_getFit(posterior,out);

Error in VBA_NLStateSpaceModel (line 364)
[posterior,out] = VBA_wrapup(posterior,options,dim,suffStat,u,y,it);

Error in demo_Qlearning (line 77)
[posterior, out] = VBA_NLStateSpaceModel(y, u, f_fname, g_fname, dim, options);

Any hints on how to solve this?

Many thanks,
Florian

It seems that some VBA functions have the same name as in SPM and interfere with SPM proper functionning. Everytime I was trying to launch my fMRI first level in SPM it wouldn't start mentionning some bug in spm_hrf.m and I struggled some time until finding where the bug came from but solved it by removing VBA from the path.

In the long-term I guess it would be good to use different names to avoid this but in the short term I'd recommend to avoid using both at the same time as they may interfere with each other.

VBA_BMA does not catch the case when some fields are missing in the posterior structure (eg a_sigma).

From @GoogleCodeExporter on October 7, 2015 9:29

This thread is a warning on the use of exceedance probabilities (EPs).

By default, a RFX-BMS will try to compute exceedance probabilities using 
compiled C code for sampling on a Dirichlet distribution (spm_gamrnd.c). If the 
compilation has not been done, this attempt will fail and the toolbox will 
resort to a normal approximation of EPs.

The quality of this approximation increases with the ratio of number of 
subjects divided by the number of models. When this ratio gets smaller than 
one, the approximation can become slightly overconfident.

Note that the toolbox prompts the user with the following message: 'Warning: 
exceedance probabilities are approximated!'.


In such cases, users would be well advised to re-compile the code, so that 
EPs are correctly estimated. In some instances, it suffices to replace the 
files 'spm_gamrnd*.*' in the directory '/stats&plots/spm_code' with the 
attached files. If this does not work (if the warning message still appears), 
please follow this thread and we will find a solution.

Original issue reported on code.google.com by [email protected] on 28 Jan 2014 at 1:53

Attachments:

• spm_gamrnd.c
• spm_gamrnd.mexmaci64
• spm_gamrnd.mexw64
• spm_gamrnd.mexa64
• spm_gamrnd.mexglx
• spm_gamrnd.mexmaci
• spm_gamrnd.mexw32

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

According to the wiki, one way to get a design allowing to better distentangle two models is to minimize the Laplace-Chernoff risk. However, it seems that VBA efficiency is "minus the Chernoff risk ". Am I right in telling that one has to maximise this output?

If it be so, could we write it explicitely to avoid any confusion?

Bastien

From @GoogleCodeExporter on October 7, 2015 9:29

What steps will reproduce the problem?
1. Fit Extended RL model (6 sessions)
2. VBA_ReDisplay_sessions(posterior, out, {[1 2]}
3. The graphical output window appears with the Summary tab looking fine, no 
error
4. When clicking on any of the other tabs ('VB inv, se...', ..., 'priors, 
ses...'), the tab is blank and an Error is thrown : 

??? Reference to non-existent field 'i_abs'.

Error in ==> VBA_ReDisplay_sessions>myPriors at 170
i_tab = ud.i_abs;

Error in ==> spm_uitab>doChoose at 195
            feval(ud.callback);

??? Error while evaluating uicontrol Callback

Toolbox updated on 2013-10-11
Matlab 7.10.0 (R2010a) 32 bits
Debian 7.1 (wheezy) 32 bits

Original issue reported on code.google.com by [email protected] on 13 Oct 2013 at 6:02

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

Hi Jean,

I have noticed that in the simulation scheme (simulateNLSS), the evolution noise enters the system when t increases 1, no matter how small the integration step deltat is. Is it the same for the inference scheme as well? Is the inference for the evolution noise adjueted for different step size deltat, (as the variance of the evolution noise should be linearly increasing with time), or the toolbox assume the step size is one when it is dealing with the evolution noise?

Thanks very much! I am looking forward to your reply!

Best,
Yan

Hi Jean,

I want to explain choices between two options and a bias existing in those choices. To do that, I have a classical softmax with the two values of the options :
𝑃(A )=1/(𝑒^((−(𝑉𝐴−𝑉𝐵))/𝛽)⁡)
In order to find how the bias toward the option A is implemented, I defined a full model with 3 more parameters:
𝑃(A )=1/(𝑒^((−(γa* 𝑉𝐴−γb* 𝑉𝐵+𝒅))/𝛽)⁡) with these hypothesis :

• the observed bias is explained by an additive value to the difference between the two options (d) (prior mean =0)
• the observed bias is explained by a multiplicative effect of one or(/and) the other option (γa and γb) (prior mean=1)
• the observed bias is explained by both a multiplicative effect and an additive effect.

I inverted the 8 possible models with your toolbox and the winning model of the comparison with VBA_groupBMC is the model with only the parameter d, with a posterior probability of 1.

When I make the comparison by computing the log-evidence of the reduced models using Savage Dickey ratios (VBA_SavageDickey), I found a totally different result : now it's the model with only γa and γb allowed to vary that win the comparison, again with a posterior probability of 1.

I can't figure out why those two kinds of inversion give me so different results. Do you have an idea ? I this not the way to test those hypothesis ? What result should I take into account ?

If you need more information, let me know,

Best,

This is to warn readers that the 'quick demo' WIKI page has a typo: the Volterra decomposition is not shown (as it should be); Rather, the comparison of simulated and estimated parameters is shown (twice).

The computation og the error term dy2 in VBA_Iphi does not exclude the points that should be excluded according to isYout. This causes bugs when NaNs are present in the data and, more problematically, can silently bias the inference.
VBA_Iphi_extended does not have the problem.

VBA_displayGroupBMC shows the exceedance probability of all models in the lower right panel.
Shouldn't it rather display the protected exceedance probability, when available?

Dear all,

I have a problem concerning the input of the VBAgroupBMC function.

the note in this function states that “L: Kxn array of log-model evidences”.

therefore, it seems that I should just put the raw glm.LogLikelihood (from fitglm in matlab). But I was instructed by a friend that I should reverse the sign of the glm.LogLikelihood before putting them in the VBAgroupBMC function. She linked to a discussion about reversing the sign of AIC (https://muut.com/i/vba-toolbox/questions:vba-groupbmc-script).

but the corresponding BMC result was inconsistent with the result that just compared group-sum BIC when I submitted the reversed glm.LogLikelihood.
actually, these two results were in the opposite direction. thus I began to reflect that maybe I shouldn't reverse the glm.LogLikelihood.

my questions is, which one should I put into the VBA_groupBMC function, the raw glm.LogLikelihood or the reversed one?

Best,

Xiaoyu

Hi,
I'm using VBA_MFX to estimate choice model parameters. Now I tapped into the following issue: I ran multiple simulations and VBA_MFX fits in a row (~ 360 times to test different models and a range of parameters). I noticed that the VBA_MFX inversion is getting slower and slower (starting from a few seconds up to >60min per run). If I restart MATLAB, the problem is gone. 'clear all' or 'close all' do not solve the problem.

Best,
Antonius

Dear all,

I found that, when using only Gaussian source in prepare_fullDCM, line 99 throws an error, as (1:sum([sources.type]~=0)) returns an empty 1x0 matrix:

inG.indr = inG.ind2(end) + (1:sum([sources.type]~=0));

Additionally, there will be no unique offset parameters for other Gaussian observational sources.

Best,
Simon

From @GoogleCodeExporter on October 7, 2015 9:29

What steps will reproduce the problem?
1. Fit Extended RL model (6 sessions)
2. Try VBA_ReDisplay_sessions(posterior, out, {[1 2]}
3. Error :
>> VBA_ReDisplay_sessions(posterior, out, {[1 2]})
Deriving diagnostics ...??? Input argument "options" is undefined.

Error in ==> VBA_LMEH0 at 15
if options.extended

Error in ==> VBA_ReDisplay_sessions>getDiagnostics at 740
[LLH0] = VBA_LMEH0(y);

Error in ==> VBA_ReDisplay_sessions at 83
ud.diagnostics = getDiagnostics(posterior,out);

-> Proposed patch : 
VBA_ReDisplay_sessions.m
740c740
< [LLH0] = VBA_LMEH0(y, out.options);

---
> [LLH0] = VBA_LMEH0(y);

=> Makes the graphical output appear with summary tab (but other tabs won't 
show - see separate issue)

Toolbox updated on 2013-10-11
Matlab 7.10.0 (R2010a) 32bits
Debian 7.1 32bits

Original issue reported on code.google.com by [email protected] on 13 Oct 2013 at 5:55

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

Some functions of the toolbox are not compatible with multisource and will fail, bug, or display inconsistent results if called on multisource models, often because if only relies on the options.binomial flag instead of using the options.source structure:

• VBA_designEfficiency
• VBA_EKF
• VBA_getDefault
• VBA_getMCMC_predictive_density_fb
• VBA_GN
• VBA_Initialize
• VBA_MFX
• VBA_onlineWrapper
• VBA_ReDisplay
• VBA_summary
• VBA_summaryMFX
• VBA_VolterraKernels

I started a branch fix-binomial to tackle those issues.

Moreover, some redundancy in the code could be avoided, as multisource coded can sometimes handle subcases that are currently treated by separate function (VBA_IPhi_XXX in particular).
Also, some features are currently incompatible with multisource although there is no principled reason for this except for the time necessary for implementing them in multisource, namely the one associated with the flag options.UNL and options.nmog (not mentioning stochastic inversion).

Dear all,

I found a small issue in VBA_bDCM_lesion in the two simulateNLSS definitions (line 18 and 37, or more specifically 25 and 44):

[yp,~,~,~,er] = simulateNLSS(...
        out.dim.n_t,...
        out.options.f_fname,...
        out.options.g_fname,...
        posterior.muTheta,...
        posterior.muPhi,...
        out.u,...
        Inf, Inf,out.options,...
        posterior.muX0);
   
        rnormal.y=yp-er;

The second Inf causes an error as simulateNLSS expects a sigma for each Gaussian source. I am really unsure about the parameters that should be fed into the function, but something like posterior.a_sigma * Inf could be an easy fix.

Best,
Simon

Hi Lionel,

The calculation of AIC and BIC values (at least when fitting multinomial data) does not adjust to the option updateX0, so that hidden states are still counted as free parameters even when they cannot be updated.

It is a very minor problem though as they can be recalculated post-hoc using LL.

Romain

I have been rerunning some inversion with the toolbox. The posteriors, fits data are the same, only the F values (initial and final) are different. They differ quite a lot on an individual subject (~4), and at the level of the group it completely changes the results. What have been changed in the computaiton of F?

Just a note but it might be worth mentionning it in the Wiki section on multisession inversion:

If you invert multisession learning models and want to skip the evaluation at state x0 (since you don't have a feedback for the first choice) you should set up options.skipfs mainly as a vector with ones at the beginning of each session. Puting options.skipfs = 1 will just skip the first observation for the first session and other estimates would be biased.

Dear all,

I realized for bDCM and DCM that Volterra Kernels are not derived for mixed observations (i.e. y matrices containing NaNs).
Using VBA_nanvar in VBA_getVolterraKernels could remove this issue (see pull request #71). But I do not know if that leads to an invalid derivation of kernels...

Best,
Simon

Hi,
VBA_MFX uses only one options.dim.n_t for all subjects. This is problematic when subjects have a different length of the experiment (e.g. due to missed trials etc.).

Following up the discussion started in issue #28

There was a major modification of the graphical engine in the 2014b release. Thoses modifications affect in particular the way "panes" are handled in the hierarchy of graphical obejcts included in the figure. Critically, if you using the wrong system leads to the plots beeing drawn beneath the pane, hence the "white page" effect. In fact the plots are there, hidden behing the white pane...

I applied a patch using a subfunction "getPanes" designed to get the handle to the selected pane irrespective of the Matlab version. Incidently, the patch does not appear in VBA_initDisplay but only in VBA_initDisplay_extended (that's why I hate duplicated code...). Therefore: I have a problem with the stochastic demo but not you, as they rely on VBA_initDisplay which works with the old panes system. You have problem with multisource demo but not I as they rely on VBA_initDisplay_extended which works with the new panes. having said that, I realize that my getPanes trick does not work for Matlab 2013a (but is fine for 2014a).

Line 358 in VBA_updateDisplay generates an error whenever options.isYout contains non-zero values (for obvious reasons).
Commenting the line solves the problem.

Hello,
after synchronizing my VBA version with the latest version on GitHub, there is a bug when using the multisource option. It seems to me that the "VBA_check" subfunction doesn't tolerate anymore to have different hyperpriors for different sources. The "demo_multisource" script does'nt work neither.
I am working Under MATLAB 2017a but the bug is also present in MATLAB 2014b.

The error message is the following:
Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in VBA_check (line 137)
options.isYout(options.sources(gsi(i)).out,t) = ~diQ;

Error in simulateNLSS (line 87)
[options,u,dim] = VBA_check(zeros(dim.p,dim.n_t),u,f_fname,g_fname,dim,options);

Error in demo_multisource (line 45)
[y,x,x0,eta,e] = simulateNLSS(T,[],g_fname,[],phi,u,Inf,sigma,options,[]);

NicoB

Hi VBA developer,
I guess this is a minor glitch, but I found that the model convergence flag, out.CV specified line 337-344, is not saved in the final 'out' object because it is not specified in the 'VBA_wrapup.m' used in line 364.
Can you verify this?

I am using the version of VBA specified below, which I get from 'out.toolbox' :
version: 'master/aa4657346cd0bda23e257bb576aa73a4c1598d52'
git: 1

Thanks,
Jungmin

Dear all,

I was testing in Behavioural DCM what happens if only the Cr matrix is defined. This leads to the following error:

Deriving prior's sufficient statistics ... 0.00 %
Error: could not initialize VB scheme : check function "easykron" at line 96

The error happens in f_dcm_extension
and happening in line 66.
dfdp(indhC,:) = inF.dhC'*easykron(ut,n);

Easykron, however, needs 3 input parameters:

function B=easykron(X,nx,n)
B = zeros(n*nx,n);
  idx = 1:n:(n*nx+eps);
  for k = 1:n  
    B(idx,k) = X;
    idx = idx+1;
  end

Unfortunately I am really unsure about how to fix it, as (to be frank) I have only a slight understanding about what is happening there.
Just by looking at the variables shape, I think something like: easykron(ut, nu, nr) could be a fix.

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 [email protected] on 26 Nov 2013 at 1:42

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

A set of functions in the root directory are undocumentemented:

• VBA_dcmMatrices (7219f27)
• VBA_get_dL (ead1bc3)
• VBA_get_tracker (a866b85)
• VBA_multisession_expand (4541e05)
• VBA_multisession_factor (4541e05)
• VBA_optimPriors (already documented)
• VBA_setup (ef43b33)
• VBA_Shapley
• VBA_susceptibility
• VBA_version (a866b85)
• factorial_struct (355d374)
• getHyperpriors (now VBA_guessHyperpriors)
• getStateParamInput
• setInputs
• setPriors
• check_struct (81d7c9d)

Just cloned the toolbox, added to my MATLAB path, and tried demo_FHN.m. What am I missing?

Main VB inversion...
Undefined function 'plotUncertainTimeSeries' for input arguments of type 'double'.

Error in VBA_updateDisplay (line 285)
            plotUncertainTimeSeries(gx,vy,dTime,display.ha(1));

Error in VBA_updateDisplay (line 34)
        VBA_updateDisplay(posterior,suffStat,options,y,it,'precisions')

Error in VBA_NLStateSpaceModel (line 255)
VBA_updateDisplay(posterior,suffStat,options,y,0,'precisions')

Error in demo_FHN (line 67)
[posterior,out] = VBA_NLStateSpaceModel(y,u,f_fname,g_fname,dim,options);

From @GoogleCodeExporter on October 7, 2015 9:29

What steps will reproduce the problem?
1. SVN checkout
2. list directories
3.

What is the expected output? What do you see instead?
Well, it should work.


What version of the product are you using? On what operating system?
I'm using version 0.0


Please provide any additional information below.
No.

Original issue reported on code.google.com by [email protected] on 26 Mar 2012 at 3:49

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