Dear Spiro,

We found a bug in the spm_BMS_gibbs.m function – the update was made in revision 4376.

As this fits between your version numbers I think it could be this.

Have a go at using your older version of SPM but with spm_BMS_gibbs.m from the newer version.

If this turns out not to be the case then can you send me your model evidence table (a matrix of subjects by models with entries being the log evidence for the relevant entry) and I’ll have a more detailed look.

I’m sorry for the inconvenience caused here.

All the Best, Will.

Dear Matthias,

Here's the updated function.

Best, Will.

________________________________________

From: Matthias Schurz [[log in to unmask]]

Sent: 02 December 2011 13:20

To: Penny, William

Subject: spm_BMS_gibbs.m - bugfix request

Dear Dr. Penny,

I read in the SPM mailing-list that there is a new fixed version of spm_BMS_gibbs.m function. could you pleas send me a copy of that? i am using bayesian model comparison for fMRI based DCM.

Thank you very much,

Matthias Schurz

--

________________________________________________

Mag.Dr. Matthias Schurz

Department of Psychology

University of Salzburg, Austria.

Tel: +43-662-8044-5142 Fax: +43-662-8044-5126

e-mail: [log in to unmask]

http://sites.google.com/site/matthiasschurz/

Am 23.06.2011 17:36, schrieb Will Penny:

> Dear all,

>

> We have recently found a bug in the spm_BMS_gibbs.m function (see below).

>

> The new fixed function will appear in the next set of bug fix updates

> - in the mean time if you'd like the new version please email me or

> Guillaume.

>

> Best,

>

> Will.

>

>

> -------- Original Message --------

>

> Dear Jan,

>

> Many thanks for that.

>

> Yes - it was a bug !

>

> I have now updated the function spm_BMS_gibbs.m - see attached.

> I have added new lines 52-58 - which ensure all log evidence

> differences are within machine range.

>

> This is the approach used in a similar function spm_BMS.m.

>

> Your colleagues' results are now changed quite a bit.

>

> 7 out of 24 subjects were effectively being ignored due to this

> overflow issue. 5 of these 7 favour model 4, and 2 favour model 2 (I

> looked in the variable g_post to find this out), with the overall

> result that now model 4 also seems a good one.

>

> Thanks for your help in tracking this down.

>

> On a *separate issue* note that if you have a different model

> comparison set (eg 5 models instead of 7) then the ordering of models

> can be different. This is a known phenomenon, a bit like the

> Alternative Vote

> (AV) system. For a concrete example have a look at section `Comparison

> Set' in the Comparing Families of DCMs paper.

>

> Best wishes,

>

> Will.

>

> On 23/06/2011 14:32, Jan Drugowitsch wrote:

>> Dear Will,

>>

>>> Can you email me the data for which this bug appears ?

>>>

>>> eg. the log model evidence matrix (lme) that is passed to

>>> spm_BMS_gibbs.

>>

>> attached is a file with 2 lme matrices, lme1 and lme2.

>>

>> lme1: 24 subjects, 7 models

>> lme2: 24 subjects, 5 models (dropping the last 2 of lme1)

>>

>> For the lme1 we get a model ordering 1>  2>  3>  ..., and for lme2 we

>> get an ordering 1>  3>  2>  ...

>>

>> In both cases, u will be infinite for at least one model/subject

>> combination due to a too large spread in lme(i, :). For this reason,

>> removing the mean is not sufficient to avoid exp() to return infinity.

>>

>> Best regards,

>> Jan

>>

>>

>>> On 23/06/2011 11:28, Guillaume Flandin wrote:

>>>>

>>>> Dear Will,

>>>>

>>>> can you have a look into this?

>>>>

>>>> Many thanks,

>>>> Guillaume.

>>>>

>>>>

>>>> -------- Original Message --------

>>>> Subject: Potential bug in DCM code

>>>> Date: Thu, 23 Jun 2011 12:21:16 +0200

>>>> From: Jan Drugowitsch<[log in to unmask]>

>>>> To: [log in to unmask]

>>>>

>>>> Hi all,

>>>>

>>>> My colleague recently tried to fit his fMRI data using DCM (RFX, no

>>>> family comparison) and ran into a weird effect that made the

>>>> relative ordering of the Bayesian model posterior change, depending

>>>> on which models where included/excluded. This obviously can't be

>>>> correct, and I traced the problem to a potential bug in

>>>> spm_BMS_gibbs.m (SPMv8, all available bug fixes applied).

>>>>

>>>> The problem seems to be in line 61:

>>>>

>>>> u=exp(lme(i,:)+log(r))+eps

>>>>

>>>> Here, if either entry of lme(i,:) takes a too large value, then the

>>>> exp() will be infinite, causing g in line 62 to be NaN.

>>>> Consecutively, all results for the corresponding subject will be

>>>> ignored without notifying the user.

>>>>

>>>> Line 49 apparently tries to avoid this problem by removing the mean

>>>> lme per user:

>>>>

>>>> lme=lme-mean(lme,2)*ones(1,Nk);

>>>>

>>>> However, this will still fail if the spread of lme's is too large.

>>>>

>>>>

>>>> I propose two ways to fix this:

>>>>

>>>> 1) remove max rather than mean:

>>>>

>>>> Line 49 could be changed to

>>>>

>>>> lme=lme-max(lme,[],2)*ones(1,Nk);

>>>>

>>>> This would guarantee that each row of lme is at most 0, avoiding

>>>> the

>>>> exp() problem at least with respect to lme. I don't know if other

>>>> parts of the SPM code use the same mean-removal trick. If they do,

>>>> then they should also be changed to max, which is safer.

>>>>

>>>> 2) remove max before taking each exp:

>>>>

>>>> Solution 1 only concerns lme, but not the second term, "log(r)" in

>>>> the exp(). So, if r is very small, log(r) will be a large negative

>>>> number, possibly causing the exp() to be inaccurate. In order to

>>>> avoid this, one could replace line 61 by the two lines

>>>>

>>>> u = lme(i,:)+log(r);

>>>> u = exp(u - max(u)) + eps;

>>>>

>>>> This is probably the safest approach. I don't think that the

>>>> additional eps is still required.

>>>>

>>>>

>>>> I hope that no other researches stumbled over the above without

>>>> noticing, which would cause their results to be wrong.

>>>>

>>>> Best regards,

>>>> Jan

>>>>

>>>>

>>>>

>>>>

>>>

From: SPM (Statistical Parametric Mapping) [mailto:[log in to unmask]] On Behalf Of Spiro Pantazatos
Sent: 02 May 2012 06:04
To: [log in to unmask]
Subject: [SPM] BMS: same data, different SPM8 versions and VERY different results

Hello SPMers,
I would be extremely appreciative if someone could help me figure out why running BMS (random effects group BMS using exceedance probabilities) on the exact same data (15 subjects: 21 models each estimated using SPM8 r4667) in two different versions of SPM8 (Version 1: SPM8 r4010; Version 2: SPM8 r4667) leads to very different results (please see attached graphs).

I ran a diff analysis between the two versions of spm_run_bms_dcm (v3995 vs. v4229) and didn't see any differences that are likely to be the issue. When I updated SPM8 from r4010 to r4667 on my old computer I was able to replicate the same results I was getting with r4667 on my other computer, and it's unlikely that any changes I made to the SPM8 code myself under the r4010 computer is the culprit since I was seeing the r4010 results on two separate installations of SPM8.

If anyone has seen this before and has any idea what could be causing this I would be extremely grateful.
Much appreciation in advance,
spiro