spm_vb_robust.m 1.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657
  1. function [slice] = spm_vb_robust (Y,slice)
  2. % Robust GLM modelling in a slice of fMRI
  3. % FORMAT [slice] = spm_vb_robust (Y,slice)
  4. %
  5. % Y - [T x N] time series with T time points, N voxels
  6. %
  7. % slice - data structure containing fields described in spm_vb_glmar.m
  8. %
  9. % Requires the 'mixture' toolbox: fullfile(spm('Dir'),'toolbox','mixture')
  10. %__________________________________________________________________________
  11. %
  12. % Reference:
  13. % W.D. Penny, J. Kilner and F. Blankenburg. Robust Bayesian General Linear
  14. % Models. NeuroImage, 36(3):661-671, 2007.
  15. %__________________________________________________________________________
  16. % Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging
  17. % Will Penny
  18. % $Id: spm_vb_robust.m 6079 2014-06-30 18:25:37Z spm $
  19. [T, Nv] = size(Y);
  20. X = slice.X;
  21. k = slice.k;
  22. for i=1:Nv,
  23. fprintf('Analysing voxel %d out of %d\n',i,Nv);
  24. yv=Y(:,i);
  25. if sum(diff(yv).^2) > 0
  26. rglm1 = spm_rglm(yv,X,1);
  27. rglm2 = spm_rglm(yv,X,2);
  28. logbf21 = rglm2.fm - rglm1.fm;
  29. w = rglm2.posts.w_mean;
  30. lambda = 1 / rglm1.variances;
  31. w_dev = sqrt(diag(rglm2.posts.w_cov));
  32. w_cov = rglm2.posts.w_cov;
  33. gamma = rglm2.posts.gamma(2,:)';
  34. else
  35. w = zeros(k,1);
  36. w_dev = w;
  37. w_cov = eye(k,k);
  38. lambda = 0;
  39. logbf21 = 0;
  40. gamma = zeros(T,1);
  41. end
  42. % Save to slice structure
  43. slice.wk_mean(:,i) = w;
  44. slice.w_dev(:,i) = w_dev;
  45. slice.w_cov{i} = w_cov;
  46. slice.mean_lambda(i) = lambda;
  47. slice.F(i) = logbf21;
  48. slice.gamma(:,i) = gamma;
  49. slice.b(:,i) = repmat(rglm2.mean_alpha, slice.k, 1);
  50. end
  51. slice.N=Nv;