nt_dprime_old

[d,err_rate,perm]=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

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

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