spm_dcm_compare.m 3.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116
  1. function spm_dcm_compare(P)
  2. % Compare two or more estimated models
  3. % FORMAT spm_dcm_compare(P)
  4. %
  5. % P - a char or cell array of DCM filenames
  6. %__________________________________________________________________________
  7. % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
  8. % Klaas Enno Stephan
  9. % $Id: spm_dcm_compare.m 3363 2009-09-04 15:11:19Z christophe $
  10. % Get DCM filenames
  11. %--------------------------------------------------------------------------
  12. if ~nargin
  13. P = spm_select([2 inf],'^DCM.*\.mat$','select DCM*.mat files');
  14. end
  15. if ~iscell(P), P = cellstr(P); end
  16. num_models = numel(P);
  17. % Load anc check all models and compute their evidence
  18. %--------------------------------------------------------------------------
  19. name = {};
  20. for model_index=1:num_models
  21. try
  22. load(P{model_index},'DCM','-mat');
  23. catch
  24. error('Cannot load DCM file "%s".',P{model_index});
  25. end
  26. % Check that none of the models is an averaged DCM
  27. %----------------------------------------------------------------------
  28. if isfield(DCM,'averaged') && DCM.averaged
  29. str = {...
  30. ['Model ' P{model_index} ' is an averaged DCM. '],...
  31. 'Please note that model comparison is not valid for averaged DCMs. ',...
  32. 'Procedure aborted.'};
  33. spm('alert*',str,'DCM <Compare>',spm('Cmdline'));
  34. return
  35. end
  36. % Check that all models refer to the same set of VOIs
  37. %----------------------------------------------------------------------
  38. if model_index == 1
  39. VOIs = {DCM.xY.name};
  40. else
  41. if ~isequal(VOIs,{DCM.xY.name})
  42. str = {...
  43. 'Selected models contain different sets of VOIs!',...
  44. 'Please note that model comparison is only valid for models with identical time series (VOIs).',...
  45. 'Procedure aborted.'};
  46. spm('alert*',str,'DCM <Compare>',spm('Cmdline'));
  47. return
  48. end
  49. end
  50. % Compute Model Evidence
  51. %----------------------------------------------------------------------
  52. [p,filename] = fileparts(P{model_index});
  53. name{end + 1} = filename;
  54. evidence(model_index) = DCM.F;
  55. end
  56. % compute conditional probability of DCMs under flat priors.
  57. %--------------------------------------------------------------------------
  58. F = evidence - min(evidence);
  59. i = F < (max(F) - 32);
  60. P = F;
  61. P(i) = max(F) - 32;
  62. P = P - min(P);
  63. P = exp(P);
  64. P = P/sum(P);
  65. % display results
  66. %--------------------------------------------------------------------------
  67. Fgraph = spm_figure('GetWin','Graphics'); figure(Fgraph); clf
  68. subplot(2,1,1)
  69. n = length(name);
  70. barh(1:n,F), hold on
  71. plot([1 1]*(max(F) - 3),[0 n],'LineWidth',4,'Color','r'), hold off
  72. set(gca,'YTick',1:n)
  73. set(gca,'YTickLabel',name)
  74. xlabel('log-evidence (relative)')
  75. title('Bayesian model comparison','FontSize',16)
  76. axis square
  77. grid on
  78. subplot(2,1,2)
  79. barh(1:n,P)
  80. set(gca,'YTick',1:n)
  81. set(gca,'YTickLabel',name)
  82. title({'conditional model probability';'under uniform model priors'})
  83. xlabel('posterior probability')
  84. axis square
  85. grid on
  86. % Output model comparison details to MATLAB command window
  87. %--------------------------------------------------------------------------
  88. assignin('base','DCM_cond_prob',P)
  89. assignin('base','DCM_log_ev',F)
  90. fprintf('\n\n DCM Model Comparison \n\n Relative Log Model Evidences:')
  91. for i = 1:length(F)
  92. fprintf( '\n %s: %-6.2f ',name{i},F(i));
  93. end
  94. fprintf('\n\n Conditional Model Probabilities under flat priors:')
  95. for i = 1:length(F)
  96. fprintf( '\n %s: %-6.2f',name{i},P(i));
  97. end
  98. fprintf('\n\n These have been stored in the workspace as ''DCM_log_ev '' and ''DCM_cond_prob''\n')