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
>>>
>>>
>>>
>>>
>>
>> --
>> William D. Penny
>> Wellcome Trust Centre for Neuroimaging
>> University College London
>> 12 Queen Square
>> London WC1N 3BG
>>
>> Tel: 020 7833 7475
>> FAX: 020 7813 1420
>> Email: [log in to unmask]
>> URL: http://www.fil.ion.ucl.ac.uk/~wpenny/
>>
>>
--
William D. Penny
Wellcome Trust Centre for Neuroimaging
University College London
12 Queen Square
London WC1N 3BG
Tel: 020 7833 7475
FAX: 020 7813 1420
Email: [log in to unmask]
URL: http://www.fil.ion.ucl.ac.uk/~wpenny/
|