Hi Rik,
It will only effect the RFX option.
With FFX the evidences are simply added up - there is no need for the
Gibbs sampling to estimate model frequencies.
Best, Will.
On 23/06/2011 17:56, Rik Henson wrote:
>
> Hi Will -
>
> In this script below, would this Gibbs error only effect 'RFX' option, or 'FFX' option too (or all options)?
>
> clear matlabbatch
> matlabbatch{1}.spm.stats.bms.bms_dcm.dir = cellstr(data_path);
> for ses=1:2
> dcmfile = {};
> for m=1:length(Model)
> dcmfile{m} = fullfile(data_path,DCMname{m}{ses});
> end
> matlabbatch{1}.spm.stats.bms.bms_dcm.sess_dcm{1}(ses).mod_dcm = cellstr(strvcat(dcmfile));
> end
> matlabbatch{1}.spm.stats.bms.bms_dcm.model_sp = {''};
> matlabbatch{1}.spm.stats.bms.bms_dcm.load_f = {''};
> matlabbatch{1}.spm.stats.bms.bms_dcm.method = 'FFX';
> %matlabbatch{1}.spm.stats.bms.bms_dcm.family_level(1).family.family_models = [];
> %matlabbatch{1}.spm.stats.bms.bms_dcm.family_level(1).family.family_names = [];
> %matlabbatch{1}.spm.stats.bms.bms_dcm.family_level(1).bma.bma_no = 0;
> matlabbatch{1}.spm.stats.bms.bms_dcm.verify_id = 0;
> spm_jobman('run',matlabbatch);
>
>
>
> -------------------------------------------------------
> Dr Richard Henson
> Assistant Director, Neuroimaging
> MRC Cognition& Brain Sciences Unit
> 15 Chaucer Road
> Cambridge, CB2 7EF, UK
>
> Office: +44 (0)1223 355 294 x522
> Mob: +44 (0)794 1377 345
> Fax: +44 (0)1223 359 062
>
> http://www.mrc-cbu.cam.ac.uk/people/rik.henson/personal
> -------------------------------------------------------
>
>
>> -----Original Message-----
>> From: SPM (Statistical Parametric Mapping) [mailto:[log in to unmask]] On
>> Behalf Of Will Penny
>> Sent: 23 June 2011 16:36
>> To: [log in to unmask]
>> Subject: [SPM] Bug in spm_BMS_gibbs.m
>>
>> 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/
>
>
--
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/
|