0001 function [d,err_rate,perm_rate,area_roc,roc]=nt_dprime(x,y,jd_flag)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 NSTEPS=1000;
0016 P=0.05;
0017 NPERMUTE=1000;
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
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;
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;
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
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
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
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
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