feComputeEvidence

 Computes a series of distance metrics between two RMSE distribtions.

 Compute summary statistics on to characterize the lesion,:
 We compute:
 - The strength of evidence, namely the effect size of the lesion
 - The Kullback-Leibler Divergence
 - Jeffrey's Divergence
 - The Eath Mover's distance

 Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.

0001 function se = feComputeEvidence(rmse1,rmse2)
0002 % Computes a series of distance metrics between two RMSE distribtions.
0003 %
0004 % Compute summary statistics on to characterize the lesion,:
0005 % We compute:
0006 % - The strength of evidence, namely the effect size of the lesion
0007 % - The Kullback-Leibler Divergence
0008 % - Jeffrey's Divergence
0009 % - The Eath Mover's distance
0010 %
0011 % Copyright (2013-2014), Franco Pestilli, Stanford University, pestillifranco@gmail.com.
0012 
0013 % Prepare the distribution of errors and the histograms describing the
0014 % erros.
0015 se.nolesion.rmse.all  = rmse1;
0016 se.lesion.rmse.all    = rmse2;
0017 se.nolesion.rmse.mean = mean(se.nolesion.rmse.all);
0018 se.lesion.rmse.mean   = mean(se.lesion.rmse.all);
0019 
0020 % Histograms
0021 se.xrange(1) = ceil(min([se.lesion.rmse.all se.nolesion.rmse.all]));
0022 se.xrange(2) = floor(max([se.lesion.rmse.all se.nolesion.rmse.all]));
0023 se.nbins     = 60;
0024 se.bins      = linspace(se.xrange(1),se.xrange(2),se.nbins);
0025 [se.lesion.hist,   se.lesion.xhist]   = hist(se.lesion.rmse.all,  se.bins);
0026 [se.nolesion.hist, se.nolesion.xhist] = hist(se.nolesion.rmse.all,se.bins);
0027 se.lesion.hist     = se.lesion.hist ./   sum(se.lesion.hist);
0028 se.nolesion.hist   = se.nolesion.hist ./ sum(se.nolesion.hist);
0029 
0030 % Kullback-Leibler Divergence
0031 se.kl.name = sprintf('Kullback–Leibler divergence: http://en.wikipedia.org/wiki/Kullback-Leibler_divergence');
0032 tmp = se.nolesion.hist .* log2( (se.nolesion.hist) ./ (se.lesion.hist + eps) );
0033 se.kl.mean = nansum(tmp);clear tmp
0034 se.kl.std  = nan;
0035 
0036 % Jeffrey's divergence
0037 se.j.name = sprintf('Jeffrey''s divergence: http://en.wikipedia.org/wiki/Divergence_(statistics)');
0038 tmp = se.nolesion.hist .* log2( (se.nolesion.hist) ./ ((se.lesion.hist + se.lesion.hist + eps)./2)  ) + ...
0039       se.lesion.hist .* log2( (se.lesion.hist) ./ ((se.lesion.hist + se.lesion.hist + eps)./2) );
0040 se.j.mean = nansum(tmp); clear tmp
0041 se.j.std   = nan;
0042 
0043 % Earth Mover's Distance:
0044 % Note: This can be very slow and my require large amounts of memory for more than 1000 voxels
0045 fprintf('[%s] Computing the Earth Mover''s distance... \n',mfilename)
0046 se.em.name = sprintf('Earth Mover''s distance: http://en.wikipedia.org/wiki/Earth_mover''s_distance');
0047 
0048 try
0049 %if (exist('emd_mex.m','file') == 2) % Using Rubinov c-code fastest
0050     pairwiseDist = zeros(size(se.lesion.xhist,2));
0051     for i=1:size(se.nolesion.xhist,2)
0052         for j=1:size(se.lesion.xhist,2)
0053             pairwiseDist(i,j) = abs(se.lesion.xhist(i)-se.nolesion.xhist(j));
0054         end
0055     end
0056     tmp_em = emd_mex(se.nolesion.hist,se.lesion.hist,pairwiseDist);
0057 catch ME %else
0058     fprintf('[%s] Cannot find compiled c-code for Earth Movers Distance.\nUsing the slower and less reliable MatLab implementation.',mfilename)
0059     [~,tmp_em] = emd(se.nolesion.xhist',se.lesion.xhist',se.nolesion.hist',se.lesion.hist',@gdf);
0060 end
0061 se.em.mean = tmp_em;
0062 se.em.std  = nan;
0063 clear tmp_emp
0064 
0065 % Strenght of evidence (effect size)
0066 fprintf('[%s] Computing the Strength of Evidence... \n',mfilename)
0067 se.s.name = sprintf('strength of evidence, d-prime: http://en.wikipedia.org/wiki/Effect_size');
0068 se.s.nboots = 5000; 
0069 se.s.nmontecarlo = 5;
0070 se.s.nbins = 200;
0071 sizeunlesioned    = length(se.nolesion.rmse.all);
0072 nullDistributionW = nan(se.s.nboots,se.s.nmontecarlo);
0073 nullDistributionWO = nan(se.s.nboots,se.s.nmontecarlo);
0074 min_x = floor(mean([se.nolesion.rmse.all]) - mean([se.nolesion.rmse.all])*.05);
0075 max_x = ceil( mean([se.lesion.rmse.all])   + mean([se.lesion.rmse.all])*.05);
0076 
0077 for inm = 1:se.s.nmontecarlo
0078     fprintf('.')
0079     parfor ibt = 1:se.s.nboots
0080         nullDistributionW(ibt,inm)  = mean(randsample(se.nolesion.rmse.all,   sizeunlesioned,true));      
0081         nullDistributionWO(ibt,inm) = mean(randsample(se.lesion.rmse.all,sizeunlesioned,true));
0082     end
0083     
0084     % Distribution unlesioned
0085     [y(:,inm),xhis] = hist(nullDistributionW(:,inm),linspace(min_x,max_x,se.s.nbins));
0086     y(:,inm) = y(:,inm)./sum(y(:,inm));
0087     
0088     % Distribution lesioned
0089     [woy(:,inm),woxhis] = hist(nullDistributionWO(:,inm),linspace(min_x,max_x,se.s.nbins));
0090     woy(:,inm) = woy(:,inm)./sum(woy(:,inm));
0091 end
0092     
0093 se.s.mean = mean(diff([mean(nullDistributionW,1); ...
0094                           mean(nullDistributionWO,1)])./sqrt(sum([std(nullDistributionW,[],1);std(nullDistributionWO,[],1)].^2,1)));
0095 se.s.std  = std(diff([mean(nullDistributionW,1); ...
0096                           mean(nullDistributionWO,1)])./sqrt(sum([std(nullDistributionW,[],1);std(nullDistributionWO,[],1)].^2,1)));
0097 disp(' done.')
0098 
0099 % Extract the mean and error of the botstrapped disributions of mean errors
0100 y_m   = mean(y,2);
0101 ywo_m = mean(woy,2);
0102 y_e   = [y_m, y_m] + 2*[-std(y,[],2),std(y,[],2)];
0103 ywo_e = [ywo_m, ywo_m] + 2*[-std(woy,[],2),std(woy,[],2)];
0104 se.s.lesioned_e   = ywo_e;
0105 se.s.lesioned_m   = ywo_m;
0106 se.s.unlesioned_e = y_e;
0107 se.s.unlesioned_m = y_m;
0108 se.s.lesioned.xbins   = woxhis;
0109 se.s.unlesioned.xbins = xhis;
0110 se.s.min_x   = min_x;
0111 se.s.max_x = max_x;
0112 
0113 end