spm_mvb_bmc.m 3.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112
  1. function [F,P,MVB] = spm_mvb_bmc(mvb)
  2. % multivariate Bayesian model comparison (Baysian decoding of a contrast)
  3. % FORMAT [F,P,MVB] = spm_mvb_bmc(mvb)
  4. %
  5. % mvb : models to compare (file names)
  6. % F : F ratio relative to null
  7. % P : P-value relative to null
  8. % MVB : best model
  9. %__________________________________________________________________________
  10. % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
  11. % Karl Friston
  12. % $Id: spm_mvb_bmc.m 5219 2013-01-29 17:07:07Z spm $
  13. %-Get figure handles and set title
  14. %--------------------------------------------------------------------------
  15. Fmvb = spm_figure('GetWin','MVB');
  16. spm_clf(Fmvb);
  17. % get MVB results
  18. %--------------------------------------------------------------------------
  19. try
  20. mvb;
  21. catch
  22. mvb = spm_select(Inf,'mat','please select models',[],pwd,'MVB_*');
  23. end
  24. MVB = load(deblank(mvb(1,:)));
  25. MVB = MVB.MVB;
  26. % display
  27. %==========================================================================
  28. figure(Fmvb)
  29. % if there is more than one MVB get maximum F for each
  30. %--------------------------------------------------------------------------
  31. if size(mvb,1) > 1
  32. X = MVB.X;
  33. name = {MVB.name(5:end)};
  34. F = max(MVB.M.F);
  35. I = 1;
  36. % check target is the same
  37. %----------------------------------------------------------------------
  38. for i = 2:size(mvb,1)
  39. MVB = load(deblank(mvb(i,:)));
  40. MVB = MVB.MVB;
  41. if ~any(X - MVB.X)
  42. name{end + 1} = MVB.name(5:end);
  43. F(end + 1) = max(MVB.M.F);
  44. I(end + 1) = i;
  45. end
  46. end
  47. % add null and model comparison
  48. %----------------------------------------------------------------------
  49. name{end + 1} = 'null';
  50. F(end + 1) = MVB.M(1).F(1);
  51. F = F - min(F);
  52. P = exp(F - mean(F));
  53. P = P./sum(P);
  54. % load best (non-null) model
  55. %----------------------------------------------------------------------
  56. [p,i] = max(P(1:end - 1));
  57. MVB = load(deblank(mvb(I(i),:)));
  58. MVB = MVB.MVB;
  59. spm_mvb_display(MVB)
  60. % display model comparison
  61. %----------------------------------------------------------------------
  62. subplot(3,2,1)
  63. bar(F), hold on
  64. plot([0 length(F) + 1], [3 3],'r')
  65. plot([0 length(F) + 1],-[3 3],'r')
  66. plot([0 length(F) + 1], [5 5],'r:')
  67. plot([0 length(F) + 1],-[5 5],'r:')
  68. hold off
  69. set(gca,'XTickLabel',name);
  70. axis square
  71. grid on
  72. title({'log-evidence';'Model comparison'})
  73. subplot(3,2,2), cla
  74. text(0,1/2,name','FontSize',12,'FontWeight','Bold')
  75. text(.7,1/2,num2str(P',' %.3f'),'FontSize',12,'FontWeight','Bold')
  76. text(.9,1/2,num2str(F',' (%2.1f)'),'FontSize',12)
  77. axis square off
  78. title({'Posterior p-values (F)';'Model comparison'})
  79. else
  80. % display
  81. %----------------------------------------------------------------------
  82. spm_mvb_display(MVB)
  83. % probability relative to null
  84. %----------------------------------------------------------------------
  85. F = MVB.M.F(2:end) - MVB.M.F(1);
  86. P = exp(F);
  87. P = P./(P + 1);
  88. end
  89. %-Reset title
  90. %--------------------------------------------------------------------------
  91. spm('Pointer','Arrow')