spm_vb_beta.m 1.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445
  1. function [block] = spm_vb_beta(Y,block)
  2. % Variational Bayes for GLM-AR modelling in a block - Update beta
  3. % FORMAT [block] = spm_vb_beta(Y,block)
  4. %
  5. % Y [T x N] time series
  6. % block data structure (see spm_vb_glmar)
  7. %__________________________________________________________________________
  8. % Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging
  9. % Will Penny and Nelson Trujillo-Barreto
  10. % $Id: spm_vb_beta.m 6079 2014-06-30 18:25:37Z spm $
  11. if block.verbose
  12. disp('Updating beta');
  13. end
  14. N = block.N;
  15. p = block.p;
  16. % Convert from V_n to V_p
  17. for n=1:N
  18. for j=1:p
  19. a_cov_p(n,j) = block.a_cov{n}(j,j);
  20. end
  21. end
  22. switch block.priors.A
  23. case 'Discrete'
  24. for j=1:p
  25. for s=1:block.priors.S
  26. sj = (block.priors.voxel(s).i-1)*p+j;
  27. H = sum((block.a_mean(sj)-block.as(j,s)).^2+a_cov_p(block.priors.voxel(s).i,j));
  28. block.b_beta(j,s) = 1/(H/2+1./block.b_beta_prior(j,s));
  29. block.mean_beta(j,s) = block.c_beta(j,s)*block.b_beta(j,s);
  30. end
  31. end
  32. otherwise
  33. for j=1:p
  34. subblock_p = j:p:N*p;
  35. H = sum(spdiags(block.Da,0).*a_cov_p(:,j)) + block.a_mean(subblock_p)'*block.Da*block.a_mean(subblock_p);
  36. % Equation 16 in paper VB4
  37. block.b_beta(j) = 1./(H./2 + 1./block.b_beta_prior(j));
  38. block.mean_beta(j) = block.c_beta(j)*block.b_beta(j);
  39. end
  40. end