Home > NoiseTools > nt_dprime.m

nt_dprime

PURPOSE ^

[d,err_rate,perm,area_roc,roc]=nt_dprime(x,y,jd_flag) - calculate d' (discriminability) of two distributions

SYNOPSIS ^

function [d,err_rate,perm_rate,area_roc,roc]=nt_dprime(x,y,jd_flag)

DESCRIPTION ^

[d,err_rate,perm,area_roc,roc]=nt_dprime(x,y,jd_flag) - calculate d' (discriminability) of two distributions

  d: discriminablity index
  err_rate: error rate for linear discrimination
  perm_rate: rate for permuted data
  area_roc: area under ROC curve
  roc: ROC curve

  x, y: data (column vectors or matrices)
  jd_flag: apply JD first
 
 NoiseTools

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [d,err_rate,perm_rate,area_roc,roc]=nt_dprime(x,y,jd_flag)
0002 %[d,err_rate,perm,area_roc,roc]=nt_dprime(x,y,jd_flag) - calculate d' (discriminability) of two distributions
0003 %
0004 %  d: discriminablity index
0005 %  err_rate: error rate for linear discrimination
0006 %  perm_rate: rate for permuted data
0007 %  area_roc: area under ROC curve
0008 %  roc: ROC curve
0009 %
0010 %  x, y: data (column vectors or matrices)
0011 %  jd_flag: apply JD first
0012 %
0013 % NoiseTools
0014 
0015 NSTEPS=1000; % number of steps to find min of error rate
0016 P=0.05; % theshold value for permutation test
0017 NPERMUTE=1000; % number of trials for permutation test
0018 
0019 if nargin<3; jd_flag=[]; end
0020 if nargin<2; error('!'); end
0021 if iscell(x)
0022     xx=[]; yy=[];
0023     for iCell=1:numel(x);
0024         xx=[xx; x{iCell}];
0025         yy=[yy; y{iCell}];
0026     end
0027     x=xx; y=yy;
0028 end
0029 if size(x,2) ~= size(y,2); error('!'); end
0030 
0031 if jd_flag; 
0032     c0=nt_cov(nt_demean(x))+nt_cov(nt_demean(y));
0033     c1=nt_cov(mean(x)-mean(y));
0034     todss=nt_dss0(c0,c1);
0035     x=nt_mmat(x,todss);
0036     y=nt_mmat(y,todss);
0037 end
0038 
0039 d=abs(mean(x)-mean(y)) ./ sqrt((var(x)+var(y))/2);
0040 
0041 % make sure that y>x
0042 for iChan=1:size(x,2)
0043     if mean(x(:,iChan))>mean(y(:,iChan));
0044         x(:,iChan)=-x(:,iChan);
0045         y(:,iChan)=-y(:,iChan);
0046     end
0047 end
0048 
0049 if nargout>1; % error rate
0050     err_rate=[];
0051     for iChan=1:size(x,2)
0052         min_error=1;
0053         for thresh=linspace(mean(x(:,iChan)),mean(y(:,iChan)),NSTEPS);
0054             x2y=sum(x(:,iChan)>thresh);
0055             y2x=sum(y(:,iChan)<thresh);
0056             nErr=x2y+y2x;
0057             min_error=min(min_error,nErr/(size(x,1)+size(y,1)));
0058         end
0059         err_rate(iChan)=min_error;
0060     end
0061 end
0062 
0063 if nargout>2; %permutation test
0064     perm_rate=[];
0065     for iChan=1:size(x,2)
0066         for iPermute=1:NPERMUTE
0067             if rem(iPermute,100)==1; disp(iPermute); end
0068             % scramble between x and y:
0069             z=[x(:,iChan);y(:,iChan)]; 
0070             z=z(randperm(size(z,1)));
0071             xx=z(1:size(x,1));    
0072             yy=z(size(x,1)+1:end);
0073             min_error=1;
0074             % scan criterion for minimum error
0075             for thresh=linspace(mean(x(:,iChan)),mean(y(:,iChan)),NSTEPS);
0076                 x2y=sum(xx>thresh);
0077                 y2x=sum(yy<thresh);
0078                 nErr=x2y+y2x;
0079                 min_error=min(min_error,nErr/(size(x,1)+size(y,1)));
0080             end
0081             min_errors(iPermute)=min_error;
0082         end
0083         % find 5th percentile of distribution of error rates
0084         min_errors=sort(min_errors);
0085         min_errors=min_errors(1:round(NPERMUTE*P));
0086         perm_rate(iChan)=min_errors(end);
0087     end
0088 end
0089 
0090 if nargout>3; error('ROC not implemented yet'); end
0091 
0092 
0093 % test code
0094 if 0
0095     x=randn(10000,1);
0096     y=1+randn(10000,1);
0097     figure(1); clf
0098     t=-3:0.1:4;
0099     plot(t,hist(x,t));
0100     hold on;
0101     plot(t,hist(y,t), 'r');
0102     [d,e,p]=nt_dprime(x,y);
0103     disp(['d'': ', num2str(d)]);
0104     disp(['e'': ', num2str(e)]);
0105     disp(['p'': ', num2str(p)]);
0106 end

Generated on Mon 14-May-2018 09:45:18 by m2html © 2005