123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116 |
- function spm_dcm_compare(P)
- % Compare two or more estimated models
- % FORMAT spm_dcm_compare(P)
- %
- % P - a char or cell array of DCM filenames
- %__________________________________________________________________________
- % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
- % Klaas Enno Stephan
- % $Id: spm_dcm_compare.m 3363 2009-09-04 15:11:19Z christophe $
- % Get DCM filenames
- %--------------------------------------------------------------------------
- if ~nargin
- P = spm_select([2 inf],'^DCM.*\.mat$','select DCM*.mat files');
- end
- if ~iscell(P), P = cellstr(P); end
- num_models = numel(P);
- % Load anc check all models and compute their evidence
- %--------------------------------------------------------------------------
- name = {};
- for model_index=1:num_models
-
- try
- load(P{model_index},'DCM','-mat');
- catch
- error('Cannot load DCM file "%s".',P{model_index});
- end
- % Check that none of the models is an averaged DCM
- %----------------------------------------------------------------------
- if isfield(DCM,'averaged') && DCM.averaged
- str = {...
- ['Model ' P{model_index} ' is an averaged DCM. '],...
- 'Please note that model comparison is not valid for averaged DCMs. ',...
- 'Procedure aborted.'};
- spm('alert*',str,'DCM <Compare>',spm('Cmdline'));
- return
- end
-
- % Check that all models refer to the same set of VOIs
- %----------------------------------------------------------------------
- if model_index == 1
- VOIs = {DCM.xY.name};
- else
- if ~isequal(VOIs,{DCM.xY.name})
- str = {...
- 'Selected models contain different sets of VOIs!',...
- 'Please note that model comparison is only valid for models with identical time series (VOIs).',...
- 'Procedure aborted.'};
- spm('alert*',str,'DCM <Compare>',spm('Cmdline'));
- return
- end
- end
-
- % Compute Model Evidence
- %----------------------------------------------------------------------
- [p,filename] = fileparts(P{model_index});
- name{end + 1} = filename;
- evidence(model_index) = DCM.F;
-
- end
- % compute conditional probability of DCMs under flat priors.
- %--------------------------------------------------------------------------
- F = evidence - min(evidence);
- i = F < (max(F) - 32);
- P = F;
- P(i) = max(F) - 32;
- P = P - min(P);
- P = exp(P);
- P = P/sum(P);
-
- % display results
- %--------------------------------------------------------------------------
- Fgraph = spm_figure('GetWin','Graphics'); figure(Fgraph); clf
- subplot(2,1,1)
- n = length(name);
- barh(1:n,F), hold on
- plot([1 1]*(max(F) - 3),[0 n],'LineWidth',4,'Color','r'), hold off
- set(gca,'YTick',1:n)
- set(gca,'YTickLabel',name)
- xlabel('log-evidence (relative)')
- title('Bayesian model comparison','FontSize',16)
- axis square
- grid on
- subplot(2,1,2)
- barh(1:n,P)
- set(gca,'YTick',1:n)
- set(gca,'YTickLabel',name)
- title({'conditional model probability';'under uniform model priors'})
- xlabel('posterior probability')
- axis square
- grid on
- % Output model comparison details to MATLAB command window
- %--------------------------------------------------------------------------
- assignin('base','DCM_cond_prob',P)
- assignin('base','DCM_log_ev',F)
- fprintf('\n\n DCM Model Comparison \n\n Relative Log Model Evidences:')
- for i = 1:length(F)
- fprintf( '\n %s: %-6.2f ',name{i},F(i));
- end
- fprintf('\n\n Conditional Model Probabilities under flat priors:')
- for i = 1:length(F)
- fprintf( '\n %s: %-6.2f',name{i},P(i));
- end
- fprintf('\n\n These have been stored in the workspace as ''DCM_log_ev '' and ''DCM_cond_prob''\n')
|