12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697 |
- function [xCon,SPM]= spm_bayes2_logbf(SPM,XYZ,xCon,ic)
- % Compute and write log Bayes factor image
- % FORMAT [xCon,SPM]= spm_bayes2_logbf(SPM,XYZ,xCon,ic)
- %
- % SPM - SPM data structure
- % XYZ - voxel list
- % xCon - contrast info
- % ic - contrast number
- %__________________________________________________________________________
- % Copyright (C) 2010-2011 Wellcome Trust Centre for Neuroimaging
- % Will Penny
- % $Id: spm_bayes2_logbf.m 4615 2012-01-10 16:56:25Z will $
-
- %-Compound Contrast
- %--------------------------------------------------------------------------
- c = xCon(ic).c;
- kc = size(c,2);
- %-Get posterior beta's
- %--------------------------------------------------------------------------
- Nk = size(SPM.xX.X,2);
- for k=1:Nk
- beta(k,:) = spm_get_data(SPM.VCbeta(k),XYZ);
- end
- %-Get noise hyperparameters
- %--------------------------------------------------------------------------
- Nh=length(SPM.PPM.l);
- for jj = 1:Nh
- hyper(jj).l = spm_get_data(SPM.VHp(jj),XYZ);
- end
-
- %-Get posterior SD beta's
- %--------------------------------------------------------------------------
- Nk = size(SPM.xX.X,2);
- %-Loop over voxels
- %==========================================================================
- Nvoxels = size(XYZ,2);
- D = NaN(reshape(SPM.xVol.DIM(1:3),1,[]));
- prior_cov=c'*SPM.PPM.Cb*c;
- a=zeros(kc,1);
- spm_progress_bar('Init',100,'Estimating Bayes factor','');
- for v=1:Nvoxels
-
- % Reconstruct approximation to voxel wise covariance matrix
- Sigma_post = SPM.PPM.Cby;
- for jj = 1:Nh,
- % Taylor approximation to posterior covariance
- Sigma_post = Sigma_post + SPM.PPM.dC{jj}*(hyper(jj).l(v) - SPM.PPM.l(jj));
- end
-
- % Get posterior for this voxel
- post_mean = c'*beta(:,v);
- post_cov = c' * Sigma_post * c;
-
- en=a-post_mean;
- iC=inv(post_cov);
- logbf = - 0.5*spm_logdet(prior_cov)+0.5*spm_logdet(post_cov)+0.5*en'*iC*en;
-
- D(XYZ(1,v),XYZ(2,v),XYZ(3,v)) = logbf;
- if rem(v,100)==0
- % update progress bar every 100th voxel
- spm_progress_bar('Set',100*v/Nvoxels);
- end
-
- end
- xCon(ic).eidf=rank(post_cov);
- spm_progress_bar('Clear');
- %-Create handle
- %--------------------------------------------------------------------------
- Vhandle = struct(...
- 'fname', [sprintf('logbf_%04d',ic) spm_file_ext],...
- 'dim', SPM.xVol.DIM',...
- 'dt', [spm_type('float32') spm_platform('bigend')],...
- 'mat', SPM.xVol.M,...
- 'pinfo', [1,0,0]',...
- 'descrip',sprintf('Log Bayes factor con %d: %s',ic,xCon(ic).name));
- %-Write image
- %--------------------------------------------------------------------------
- Vhandle = spm_create_vol(Vhandle);
- Vhandle = spm_write_vol(Vhandle,D);
- xCon(ic).Vcon = Vhandle;
- fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),...
- sprintf('...written %s',spm_file(Vhandle.fname,'filename'))); %-#
|