123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112 |
- function [F,P,MVB] = spm_mvb_bmc(mvb)
- % multivariate Bayesian model comparison (Baysian decoding of a contrast)
- % FORMAT [F,P,MVB] = spm_mvb_bmc(mvb)
- %
- % mvb : models to compare (file names)
- % F : F ratio relative to null
- % P : P-value relative to null
- % MVB : best model
- %__________________________________________________________________________
- % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
- % Karl Friston
- % $Id: spm_mvb_bmc.m 5219 2013-01-29 17:07:07Z spm $
- %-Get figure handles and set title
- %--------------------------------------------------------------------------
- Fmvb = spm_figure('GetWin','MVB');
- spm_clf(Fmvb);
- % get MVB results
- %--------------------------------------------------------------------------
- try
- mvb;
- catch
- mvb = spm_select(Inf,'mat','please select models',[],pwd,'MVB_*');
- end
- MVB = load(deblank(mvb(1,:)));
- MVB = MVB.MVB;
- % display
- %==========================================================================
- figure(Fmvb)
- % if there is more than one MVB get maximum F for each
- %--------------------------------------------------------------------------
- if size(mvb,1) > 1
-
- X = MVB.X;
- name = {MVB.name(5:end)};
- F = max(MVB.M.F);
- I = 1;
-
- % check target is the same
- %----------------------------------------------------------------------
- for i = 2:size(mvb,1)
- MVB = load(deblank(mvb(i,:)));
- MVB = MVB.MVB;
- if ~any(X - MVB.X)
- name{end + 1} = MVB.name(5:end);
- F(end + 1) = max(MVB.M.F);
- I(end + 1) = i;
- end
- end
- % add null and model comparison
- %----------------------------------------------------------------------
- name{end + 1} = 'null';
- F(end + 1) = MVB.M(1).F(1);
- F = F - min(F);
- P = exp(F - mean(F));
- P = P./sum(P);
-
- % load best (non-null) model
- %----------------------------------------------------------------------
- [p,i] = max(P(1:end - 1));
- MVB = load(deblank(mvb(I(i),:)));
- MVB = MVB.MVB;
- spm_mvb_display(MVB)
-
- % display model comparison
- %----------------------------------------------------------------------
- subplot(3,2,1)
- bar(F), hold on
- plot([0 length(F) + 1], [3 3],'r')
- plot([0 length(F) + 1],-[3 3],'r')
- plot([0 length(F) + 1], [5 5],'r:')
- plot([0 length(F) + 1],-[5 5],'r:')
- hold off
- set(gca,'XTickLabel',name);
- axis square
- grid on
- title({'log-evidence';'Model comparison'})
- subplot(3,2,2), cla
- text(0,1/2,name','FontSize',12,'FontWeight','Bold')
- text(.7,1/2,num2str(P',' %.3f'),'FontSize',12,'FontWeight','Bold')
- text(.9,1/2,num2str(F',' (%2.1f)'),'FontSize',12)
- axis square off
- title({'Posterior p-values (F)';'Model comparison'})
- else
-
- % display
- %----------------------------------------------------------------------
- spm_mvb_display(MVB)
-
- % probability relative to null
- %----------------------------------------------------------------------
- F = MVB.M.F(2:end) - MVB.M.F(1);
- P = exp(F);
- P = P./(P + 1);
- end
- %-Reset title
- %--------------------------------------------------------------------------
- spm('Pointer','Arrow')
|