nt_dprime

[d,e]=nt_dprime(x,y,jd_flag) - calculate d' (discriminability) of two distributions

  d: discriminablity index
  e: error rate for linear discrimination

  x, y: data (observtions X features) - if cell arrays do xvalidation
  jd_flag: apply JD first
 
 See nt_dprime_old for earlier implementation (deprecated).
 
 NoiseTools

0001 function [d,e]=nt_dprime(x,y,jd_flag)
0002 %[d,e]=nt_dprime(x,y,jd_flag) - calculate d' (discriminability) of two distributions
0003 %
0004 %  d: discriminablity index
0005 %  e: error rate for linear discrimination
0006 %
0007 %  x, y: data (observtions X features) - if cell arrays do xvalidation
0008 %  jd_flag: apply JD first
0009 %
0010 % See nt_dprime_old for earlier implementation (deprecated).
0011 %
0012 % NoiseTools
0013 
0014 if nargin<3; jd_flag=[]; end
0015 if nargin<2; error('!'); end
0016 
0017 if ~iscell(x); x={x,x}; end % no crossvalidation but will use same code
0018 if ~iscell(y); y={y,y}; end
0019 
0020 for iCell=1:numel(x)
0021     if size(x{iCell},2) ~= size(y{iCell},2); error('!'); end
0022     xx=x{iCell}; yy=y{iCell};
0023     others=setdiff(1:numel(x),iCell);
0024     xxx=x(others); xxx=cat(1,xxx{:});
0025     yyy=y(others); yyy=cat(1,yyy{:});
0026     if jd_flag; 
0027         c0=nt_cov(nt_demean(xxx))+nt_cov(nt_demean(yyy));
0028         c1=nt_cov(mean(xxx)-mean(yyy));
0029         todss=nt_dss0(c0,c1);
0030         xxx=nt_mmat(xxx,todss); % others
0031         yyy=nt_mmat(yyy,todss);
0032         xx=nt_mmat(xx,todss); % this cell
0033         yy=nt_mmat(yy,todss);
0034     end
0035 
0036     % d-prime on left-out data
0037     d(iCell,:)=abs(mean(xx)-mean(yy)) ./ sqrt((var(xx)+var(yy))/2);
0038 
0039     % error rate on left-out data
0040     for iChan=1:size(xx,2);
0041         dxx=(xx(:,iChan)-mean(xxx(:,iChan))).^2;
0042         dxy=(xx(:,iChan)-mean(yyy(:,iChan))).^2;
0043         dyx=(yy(:,iChan)-mean(xxx(:,iChan))).^2;
0044         dyy=(yy(:,iChan)-mean(yyy(:,iChan))).^2;
0045         nMissX=numel(xx(find(dxx>dxy),iChan));
0046         nMissY=numel(yy(find(dyy>dyx),iChan));
0047         ex(iCell,iChan)=nMissX/numel(xx(:,iChan)); % xs misclassified as ys
0048         ey(iCell,iChan)=nMissY/numel(yy(:,iChan)); % ys misclassified as xs
0049     end
0050 end
0051 
0052 d=mean(d); 
0053 ex=mean(ex);
0054 ey=mean(ey);
0055 e=(ex+ey)/2;
0056 
0057 
0058 
0059 % test code
0060 if 0
0061     x=randn(10000,1);
0062     y=1+randn(10000,1);
0063     figure(1); clf
0064     t=-3:0.1:6;
0065     plot(t,hist(x,t));
0066     hold on;
0067     plot(t,hist(y,t), 'r');
0068     [d,e]=nt_dprime(x,y);
0069     disp(['d'': ', num2str(d)]);
0070     disp(['error: ', num2str(e)]);
0071 end
0072 if 0
0073     x={randn(10000,1),randn(10000,1)};
0074     y={1+randn(10000,1),1+randn(10000,1)};
0075     [d,e]=nt_dprime(x,y);
0076     disp(['d'': ', num2str(d)]);
0077     disp(['error: ', num2str(e)]);
0078 end