spm_BMS_F.m 1.3 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950
  1. function [F_samp,F_bound] = spm_BMS_F (alpha,lme,alpha0)
  2. % Compute two lower bounds on model evidence p(y|r) for group BMS
  3. %
  4. % FORMAT [F_samp,F_bound] = spm_BMS_smpl_me (alpha,lme,alpha0)
  5. %
  6. % INPUT:
  7. % alpha parameters of p(r|y)
  8. % lme array of log model evidences
  9. % rows: subjects
  10. % columns: models (1..Nk)
  11. % alpha0 priors of p(r)
  12. %
  13. % OUTPUT:
  14. % F_samp - sampling estimate of <ln p(y_n|r>
  15. % F_bound - lower bound on lower bound of <ln p(y_n|r>
  16. %
  17. % REFERENCE: See appendix in
  18. % Stephan KE, Penny WD, Daunizeau J, Moran RJ, Friston KJ
  19. % Bayesian Model Selection for Group Studies. NeuroImage (under review)
  20. %__________________________________________________________________________
  21. % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
  22. % Will Penny
  23. % $Id: spm_BMS_F.m 2507 2008-11-30 14:45:22Z klaas $
  24. alpha0 = sort(alpha0);
  25. if alpha0(1) ~= alpha0(end)
  26. error('Error in function spm_BMS_F: alpha0 should have identical values.')
  27. end
  28. alpha0 = alpha0(1);
  29. a_sum = sum(alpha);
  30. psi_sum = psi(a_sum);
  31. psi_diff = psi(alpha) - psi_sum;
  32. gm = gammaln(alpha);
  33. [s_samp,s_bound] = spm_BMS_F_smpl(alpha,lme,alpha0);
  34. K = length(alpha);
  35. F = 0;
  36. for k = 1:K,
  37. F = F - (alpha(k) - alpha0)*psi_diff(k) + gm(k);
  38. end
  39. F = F - gammaln(a_sum);
  40. F_bound = F + s_bound;
  41. F_samp = F + s_samp;
  42. return