1function [CC,Q,tsd,md]=findclassifier2(D,TRIG,cl,T,t0,SWITCH) 2% FINDCLASSIFIER2 3% is very similar to FINDCLASSIFIER1 but uses a different 4% criterion for selecting the time segment. 5% [CC,Q,TSD,MD]=findclassifier2(D,TRIG,Class,class_times,t_ref); 6% 7% D data, each row is one time point 8% TRIG trigger time points 9% Class class information 10% class_times classification times, combinations of times must be in one row 11% t_ref reference time for Class 0 (optional) 12% 13% CC contains LDA and MD classifiers 14% Q is a list of classification quality for each time of 'class_times' 15% TSD returns the LDA classification 16% MD returns the MD classification 17% 18% [CC,Q,TSD,MD]=findclassifier2(AR,find(trig>0.5)-257,~mod(1:80,2),reshape(1:14*128,16,14*8)'); 19% 20% 21% Reference(s): 22% [1] Schlögl A., Neuper C. Pfurtscheller G. 23% Estimating the mutual information of an EEG-based Brain-Computer-Interface 24% Biomedizinische Technik 47(1-2): 3-8, 2002. 25% [2] A. Schlögl, C. Keinrath, R. Scherer, G. Pfurtscheller, 26% Information transfer of an EEG-based Bran-computer interface. 27% Proceedings of the 1st International IEEE EMBS Conference on Neural Engineering, Capri, Italy, Mar 20-22, 2003 28 29 30% $Id$ 31% Copyright (C) 1999-2005 by Alois Schloegl <alois.schloegl@gmail.com> 32% This is part of the BIOSIG-toolbox http://biosig.sf.net/ 33 34 35% This program is free software; you can redistribute it and/or 36% modify it under the terms of the GNU General Public License 37% as published by the Free Software Foundation; either version 2 38% of the License, or (at your option) any later version. 39% 40% This program is distributed in the hope that it will be useful, 41% but WITHOUT ANY WARRANTY; without even the implied warranty of 42% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 43% GNU General Public License for more details. 44% 45% You should have received a copy of the GNU General Public License 46% along with this program; if not, write to the Free Software 47% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 48 49 50warning('this function is obsolete and replaced by FINDCLASSIFIER'); 51 52 53tsd=[];md=[]; 54 55if nargin<6, 56 SWITCH=0; 57end; 58if nargin>4, 59 if isempty(t0), 60 t0=logical(ones(size(T,1),1)); 61 end; 62end; 63tmp=cl;tmp(isnan(tmp))=0; 64if any(rem(tmp,1) & ~isnan(cl)), 65 fprintf(2,'Error %s: class information is not integer\n',mfilename); 66 return; 67end; 68if length(TRIG)~=length(cl); 69 fprintf(2,'number of Triggers do not match class information'); 70end; 71 72cl = cl(:); 73CL = unique(cl(~isnan(cl))); 74 75[CL,iCL] = sort(CL); 76TRIG = TRIG(:); 77if ~all(D(:,1)==1) 78 % D1=[ones(size(D,1)-1,1),diff(D)]; 79 D = [ones(size(D,1),1),D]; 80 %else 81 % D1=[ones(size(D,1)-1,1),diff(D(:,2:end))]; 82end; 83 84% add sufficient NaNs at the beginning and the end 85tmp = min(TRIG)+min(min(T))-1; 86if tmp<0, 87 TRIG = TRIG - tmp; 88 D = [repmat(nan,[-tmp,size(D,2)]);D]; 89end; 90tmp = max(TRIG)+max(max(T))-size(D,1); 91if tmp>0, 92 D = [D;repmat(nan,[tmp,size(D,2)])]; 93end; 94 95% estimate classification result for all time segments in T - without crossvalidation 96CMX = zeros([size(T,1),length(CL)*[1,1]]); 97for k = 1:size(T,1), 98 cmx = zeros(length(CL)); 99 for l = 1:length(CL), 100 t = perm(TRIG(cl==CL(l)),T(k,:)); 101 %t = t(t<=size(D,1)); 102 [C{k,l},tmp] = covm(D(t(:),:),'M'); 103 end; 104 %[Q(k),d{k}] = qcmahal({C0r,C{k,:}}); 105 [CC.QC(k),d{k}] = qcmahal({C{k,:}}); 106 %[CC.QC3(k,1:2),d3{k}] = qcloglik({C{k,:}}); 107 CC.QC2(k) = sum(sqrt(d{k}(:)))/(size(d{k},1)*(size(d{k},1)-1)); 108 lnQ(k) = mean(log(d{k}(~eye(length(d{k}))))); 109 for l = 1:length(CL), 110 t = perm(TRIG(cl==CL(l)),T(k,:)); 111 %t = t(t<=size(D,1)); 112 [tmp] = mdbc({C{k,:}},D(t(:),:)); 113 [tmp,ix] = min(tmp,[],2); 114 tmp = isnan(tmp); 115 ix(tmp) = NaN; %NC(1)+1; % not classified; any value but not 1:length(MD) 116 ix(~tmp) = CL(ix(~tmp)); 117 tmp = histo3([ix;CL(:)]); 118 cmx(tmp.X,l) = tmp.H-1; 119 120 [tmp] = gdbc({C{k,:}},D(t(:),:)); 121 [tmp,ix] = min(tmp,[],2); 122 tmp = isnan(tmp); 123 ix(tmp) = NaN; %NC(1)+1; % not classified; any value but not 1:length(MD) 124 ix(~tmp) = CL(ix(~tmp)); 125 tmp = histo3([ix;CL(:)]); 126 cmx2(iCL(tmp.X),l) = tmp.H-1; 127 end; 128 CMX(k,:,:) = cmx; 129 CC.KAPPA(k) = kappa(cmx); 130 CC.ACC(k) = sum(diag(cmx))/sum(cmx(:)); 131 CC.KAPPA2(k) = kappa(cmx2); 132 CC.ACC2(k) = sum(diag(cmx2))/sum(cmx2(:)); 133 134 t = perm(TRIG(~isnan(cl)),T(k,:)); 135 tmp = repmat(cl(~isnan(cl))',size(T,2),1); 136 [tmp,acc_g(k,:),kap_g(k,:),H] = gdbc({C{k,:}},D(t(:),:),tmp(:)); 137end; 138CC.GDBC.acc = acc_g; 139CC.GDBC.kap = kap_g; 140CC.GDBC.H = H; 141CC.GDBC.p = tmp; 142 143% identify best classification time 144if nargin>4, 145 tmp = CC.QC; 146 tmp(~t0) = 0; 147 [maxQ,CC.TI] = max(tmp); %d{K}, 148else 149 [maxQ,CC.TI] = max(CC.QC); %d{K}, 150end; 151CC.qc = [CC.QC',CC.KAPPA',CC.ACC',CC.QC2']; 152%plot([CC.QC',CC.KAPPA',CC.ACC',CC.QC2',CC.QC3',t0]);drawnow 153[maxQ(1),TI(1)] = max(CC.QC'.*t0); %d{K}, 154[maxQ(2),TI(2)] = max(CC.KAPPA'.*t0); %d{K}, 155[maxQ(3),TI(3)] = max(CC.ACC'.*t0); %d{K}, 156[maxQ(4),TI(4)] = max(CC.QC2'.*t0); %d{K}, 157%[maxQ(5),TI(5)] = max(CC.QC3(:,1).*t0); %d{K}, 158%[maxQ(6),TI(6)] = max(CC.QC3(:,2).*t0); %d{K}, 159CC.TIS = TI; 160CC.TI = TI(2); 161%if any(TI(1)~=TI), 162% fprintf(1,'FINDCLASSIFIER2: \nTIS: %i\t%i\t%i\t%i',TI); 163% fprintf(1,'\n %5.3f\t%5.3f\t%5.3f\t%5.3f\n',maxQ); 164%end; 165 166% build MD classifier 167CC.D = d{CC.TI}; 168CC.Q = CC.QC(CC.TI); 169CC.MD = {C{CC.TI,:}}; 170 171if SWITCH<0, return; end; 172 173CC.IR = mdbc({C{CC.TI,:}}); 174CC.CMX = squeeze(CMX(CC.TI,:,:)); 175CC.CMX = CMX; 176CC.tsc = T(CC.TI,[1,end]); 177 178m1=decovm(CC.MD{1}); 179m2=decovm(CC.MD{2}); 180tmp=mdbc(CC.MD,[1,m1;1,m2]); 181CC.scale=[1,1]*max(abs(tmp(:))); % element 1 182 183[maxQ,CC.lnTI] = max(lnQ); %d{K}, 184CC.DistMXln = d{CC.lnTI}; 185CC.MDln = {C{CC.lnTI,:}}; 186 187% alternative classifier using two different time segments. 188if 1, 189 [Q,d] = qcmahal({C{:}}'); 190 CC.T2.D = d; 191 [ix,iy] = find(d==max(d(:))); 192 ix=mod(ix(1)-1,size(C,1))+1; 193 iy=mod(iy(1)-1,size(C,1))+1; 194 CC.T2.TI = [ix(1),iy(1)]; 195 CC.T2.MD = {C{ix(1),1},C{iy(1),2}}; 196end; 197 198% build LDA classifier 199if length(CL)==2, 200 % LDA 201 C0 = zeros(size(C{CC.TI,1})); 202 for l=1:length(CL); 203 [M{l},sd,COV,xc,N,R2] = decovm(C{CC.TI,l}); 204 C0 = C0 + C{CC.TI,l}; 205 end; 206 [M0,sd,COV0,xc,N,R2] = decovm(C0); 207 w = COV0\(M{1}'-M{2}'); 208 w0 = M0*w; 209 %CC.LDA.b = w0; 210 %CC.LDA.w = -w; 211 CC.lda = [w0; -w]; 212 213 % MD 214 tsd = ldbc(CC.MD,D); 215end; 216md = mdbc(CC.MD,D); 217lld= llbc(CC.MD,D); 218 219% bias correction not used anymore. 220CC.BIAS.LDA=0;%mean(tsd); 221CC.BIAS.MDA=0;%mean(diff(-md,[],2)); 222CC.BIAS.GRB=0;%mean(diff(-exp(-md/2),[],2)); 223 224% define datatype 225CC.datatype = 'Classifier'; 226CC.Classes = CL; 227 228 229%% cross-validation with jackknife (leave one trial out) 230nc = max(max(T))-min(min(T))+1; 231 232JKD = repmat(nan,[nc,length(CL),length(TRIG)]); 233JKLL = JKD; 234JKD1 = repmat(nan,[nc,length(TRIG)]); 235JKD2 = repmat(nan,[nc,length(TRIG)]); 236JKD3 = repmat(nan,[nc,length(TRIG)]); 237JKD4 = repmat(nan,[nc,length(TRIG)]); 238JKD5 = repmat(nan,[nc,length(TRIG)]); 239JKD6 = repmat(nan,[nc,length(TRIG)]); 240JKLD = repmat(nan,[nc,length(TRIG)]); 241for l = find(~isnan(cl(:)'));1:length(cl); 242 c = find(cl(l)==CL); 243 t = TRIG(l)+T(CC.TI,:); 244 %t = t(t<=size(D,1)); 245 [tmp,tmpn] = covm(D(t(:),:),'M'); 246 247 cc = CC.MD; 248 cc{c} = CC.MD{c}-tmp; 249 250 %t = TRIG(l)+(1:nc); 251 %t = t(t<=size(D,1)); 252 t = TRIG(l)+(min(min(T)):max(max(T))); 253 254 d = llbc(cc,D(t,:)); 255 JKLL(:,:,l) = d; 256 [tmp,LLIX(:,l)] = max(d,[],2); 257 if length(CL)==2, 258 JKD3(:,l)=d(:,1); 259 JKD4(:,l)=d(:,2); 260 end; 261 262 d = mdbc(cc,D(t,:)); 263 JKD(:,:,l) = d; 264 [tmp,MDIX(:,l)] = min(d,[],2); 265 266 if length(CL)==2, 267 JKD1(:,l) = d(:,1); 268 JKD2(:,l) = d(:,2); 269 end; 270 271 d = gdbc(cc,D(t,:)); 272 JKGD(:,:,l) = d; 273 [tmp,GDIX(:,l)] = max(d,[],2); 274 275 d = ldbc2(cc,D(t,:)); 276 LDA2(:,:,l) = d; 277 [tmp,LD2IX(:,l)] = max(d,[],2); 278 279 d = ldbc3(cc,D(t,:)); 280 LDA3(:,:,l) = d; 281 size(LDA3), 282 [tmp,LD3IX(:,l)] = max(d,[],2); 283 284 d = ldbc4(cc,D(t,:)); 285 LDA4(:,:,l) = d; 286 [tmp,LD4IX(:,l)] = max(d,[],2); 287 288 if length(CL)==2, 289 JKD5(:,l) = d(:,1); 290 JKD6(:,l) = d(:,2); 291 292 LDA(:,l) = ldbc(cc); 293 JKLD(:,l) = D(t,:)*LDA(:,l); 294 end; 295end; 296if length(CL)==2, 297 [CC.ldaC0,NN] = covm(LDA','D0'); 298 %CC.ldaC0=CC.ldaC0./NN*min(0,sum(~isnan(CL))-1); 299 % since NN==min(0,sum(~isnan(CL))-1), no need to rescale 300end; 301 302% Concordance matrix with cross-validation 303CC.lmx= zeros([size(LLIX,1),length(CL)^2]); 304CC.LLH.I0 = zeros([size(LLIX,1),length(CL)]); 305CC.LLH.I = zeros([size(LLIX,1),1]); 306 307CC.mmx= zeros([size(MDIX,1),length(CL)^2]); 308CC.MD2.I0 = zeros([size(MDIX,1),length(CL)]); 309CC.MD2.I = zeros([size(MDIX,1),1]); 310tmp = zeros([size(MDIX,1),length(CL)]); 311for k = 1:length(CL), 312 313 jkd = squeeze(LDA4(:,k,:)); 314 o = bci3eval(jkd(:,(cl~=k)&~isnan(cl)),jkd(:,cl==k),2); 315 CC.LD4.TSD{k} = o; 316 CC.LD4.I0(:,k) = log2(2*var(jkd,[],2)./(var(jkd(:,cl==k),[],2) + var(jkd(:,(cl~=k)&~isnan(cl)),[],2)))/2; 317 [sum0,n0,ssq0]=sumskipnan(jkd(:,cl==k),2); 318 [sum1,n1,ssq1]=sumskipnan(jkd(:,(cl~=k)&~isnan(cl)),2); 319 s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 320 s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 321 s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 322 SNR = 2*s./(s0+s1); % this is SNR+1 323 CC.LD4.I1(:,k) = log2(SNR)/2; 324 325 jkd = squeeze(LDA3(:,k,:)); 326 o = bci3eval(jkd(:,(cl~=k)&~isnan(cl)),jkd(:,cl==k),2); 327 CC.LD3.TSD{k} = o; 328 CC.LD3.I0(:,k) = log2(2*var(jkd,[],2)./(var(jkd(:,cl==k),[],2) + var(jkd(:,(cl~=k)&~isnan(cl)),[],2)))/2; 329 [sum0,n0,ssq0]=sumskipnan(jkd(:,cl==k),2); 330 [sum1,n1,ssq1]=sumskipnan(jkd(:,(cl~=k)&~isnan(cl)),2); 331 s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 332 s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 333 s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 334 SNR = 2*s./(s0+s1); % this is SNR+1 335 CC.LD3.I1(:,k) = log2(SNR)/2; 336 337 jkd = squeeze(LDA2(:,k,:)); 338 o = bci3eval(jkd(:,(cl~=k)&~isnan(cl)),jkd(:,cl==k),2); 339 CC.LD2.TSD{k} = o; 340 CC.LD2.I0(:,k) = log2(2*var(jkd,[],2)./(var(jkd(:,cl==k),[],2) + var(jkd(:,(cl~=k)&~isnan(cl)),[],2)))/2; 341 [sum0,n0,ssq0]=sumskipnan(jkd(:,cl==k),2); 342 [sum1,n1,ssq1]=sumskipnan(jkd(:,(cl~=k)&~isnan(cl)),2); 343 s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 344 s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 345 s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 346 SNR = 2*s./(s0+s1); % this is SNR+1 347 CC.LD2.I1(:,k) = log2(SNR)/2; 348 349 jkd = squeeze(JKGD(:,k,:)); 350 o = bci3eval(jkd(:,(cl~=k)&~isnan(cl)),jkd(:,cl==k),2); 351 CC.MD3.TSD{k} = o; 352 CC.MD3.I0(:,k) = log2(2*var(jkd,[],2)./(var(jkd(:,cl==k),[],2) + var(jkd(:,(cl~=k)&~isnan(cl)),[],2)))/2; 353 [sum0,n0,ssq0]=sumskipnan(jkd(:,cl==k),2); 354 [sum1,n1,ssq1]=sumskipnan(jkd(:,(cl~=k)&~isnan(cl)),2); 355 s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 356 s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 357 s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 358 SNR = 2*s./(s0+s1); % this is SNR+1 359 CC.MD3.I1(:,k) = log2(SNR)/2; 360 361 jkd = squeeze(JKD(:,k,:)); 362 o = bci3eval(jkd(:,(cl~=k)&~isnan(cl)),jkd(:,cl==k),2); 363 CC.MD2.TSD{k} = o; 364 CC.MD2.I0(:,k) = log2(2*var(jkd,[],2)./(var(jkd(:,cl==k),[],2) + var(jkd(:,(cl~=k)&~isnan(cl)),[],2)))/2; 365 [sum0,n0,ssq0]=sumskipnan(jkd(:,cl==k),2); 366 [sum1,n1,ssq1]=sumskipnan(jkd(:,(cl~=k)&~isnan(cl)),2); 367 s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 368 s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 369 s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 370 SNR = 2*s./(s0+s1); % this is SNR+1 371 CC.MD2.I1(:,k) = log2(SNR)/2; 372 373 jkd = squeeze(JKD(:,k,:).^2); 374 o = bci3eval(jkd(:,(cl~=k)&~isnan(cl)),jkd(:,cl==k),2); 375 CC.MDA.TSD{k} = o; 376 CC.MDA.I0(:,k) = log2(2*var(jkd,[],2)./(var(jkd(:,cl==k),[],2) + var(jkd(:,(cl~=k)&~isnan(cl)),[],2)))/2; 377 [sum0,n0,ssq0]=sumskipnan(jkd(:,cl==k),2); 378 [sum1,n1,ssq1]=sumskipnan(jkd(:,(cl~=k)&~isnan(cl)),2); 379 s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 380 s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 381 s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 382 SNR = 2*s./(s0+s1); % this is SNR+1 383 CC.MDA.I1(:,k) = log2(SNR)/2; 384 385 jkd = exp(-squeeze(JKD(:,k,:))/2); 386 o = bci3eval(jkd(:,(cl~=k)&~isnan(cl)),jkd(:,cl==k),2); 387 CC.GRB.TSD{k} = o; 388 CC.GRB.I0(:,k) = log2(2*var(jkd,[],2)./(var(jkd(:,cl==k),[],2) + var(jkd(:,(cl~=k)&~isnan(cl)),[],2)))/2; 389 [sum0,n0,ssq0] = sumskipnan(jkd(:,cl==k),2); 390 [sum1,n1,ssq1] = sumskipnan(jkd(:,(cl~=k)&~isnan(cl)),2); 391 s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 392 s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 393 s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 394 SNR = 2*s./(s0+s1); % this is SNR+1 395 CC.GRB.I1(:,k) = log2(SNR)/2; 396 397 jkd = exp(-squeeze(JKLL(:,k,:))/2); 398 o = bci3eval(jkd(:,(cl~=k)&~isnan(cl)),jkd(:,cl==k),2); 399 CC.LLH.TSD{k} = o; 400 CC.LLH.I0(:,k) = log2(2*var(jkd,[],2)./(var(jkd(:,cl==k),[],2) + var(jkd(:,(cl~=k)&~isnan(cl)),[],2)))/2; 401 [sum0,n0,ssq0] = sumskipnan(jkd(:,cl==k),2); 402 [sum1,n1,ssq1] = sumskipnan(jkd(:,(cl~=k)&~isnan(cl)),2); 403 s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 404 s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 405 s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 406 SNR = 2*s./(s0+s1); % this is SNR+1 407 CC.LLH.I1(:,k) = log2(SNR)/2; 408 409 for l = 1:length(CL), 410 tmp(:,l) = sum(MDIX(:,cl==CL(k))==CL(l),2); 411 if CL(k) == CL(l), 412 acc = tmp(:,l); 413 end; 414 end; 415 CC.MDA.mmx(:,(1-length(CL):0)+k*length(CL)) = tmp; 416 CC.MDA.acc(:,k) = acc./sum(tmp,2); 417 418 for l = 1:length(CL), 419 tmp(:,l) = sum(GDIX(:,cl==CL(k))==CL(l),2); 420 if CL(k) == CL(l), 421 acc = tmp(:,l); 422 end; 423 end; 424 CC.MD3.mmx(:,(1-length(CL):0)+k*length(CL)) = tmp; 425 CC.MD3.acc(:,k) = acc./sum(tmp,2); 426 427 for l = 1:length(CL), 428 tmp(:,l) = sum(LD2IX(:,cl==CL(k))==CL(l),2); 429 if CL(k) == CL(l), 430 acc = tmp(:,l); 431 end; 432 end; 433 CC.LD2.mmx(:,(1-length(CL):0)+k*length(CL)) = tmp; 434 CC.LD2.acc(:,k) = acc./sum(tmp,2); 435 436 for l = 1:length(CL), 437 tmp(:,l) = sum(LD3IX(:,cl==CL(k))==CL(l),2); 438 if CL(k) == CL(l), 439 acc = tmp(:,l); 440 end; 441 end; 442 CC.LD3.mmx(:,(1-length(CL):0)+k*length(CL)) = tmp; 443 CC.LD3.acc(:,k) = acc./sum(tmp,2); 444 445 for l = 1:length(CL), 446 tmp(:,l) = sum(LD4IX(:,cl==CL(k))==CL(l),2); 447 if CL(k) == CL(l), 448 acc = tmp(:,l); 449 end; 450 end; 451 CC.LD4.mmx(:,(1-length(CL):0)+k*length(CL)) = tmp; 452 CC.LD4.acc(:,k) = acc./sum(tmp,2); 453 454 for l = 1:length(CL), 455 tmp(:,l) = sum(LLIX(:,cl==CL(k))==CL(l),2); 456 if CL(k) == CL(l), 457 acc = tmp(:,l); 458 end; 459 end; 460 CC.LLH.mmx(:,(1-length(CL):0)+k*length(CL)) = tmp; 461 CC.LLH.acc(:,k) = acc./sum(tmp,2); 462end; 463CC.MDA.I = sum(CC.MDA.I0,2); 464CC.MD2.I = sum(CC.MD2.I0,2); 465CC.MD3.I = sum(CC.MD3.I0,2); 466CC.LD2.I = sum(CC.LD2.I0,2); 467CC.GRB.I = sum(CC.GRB.I0,2); 468CC.LLH.I = sum(CC.LLH.I0,2); 469 470CC.MDA.CMX00 = reshape(sum(CC.MDA.mmx(T(CC.TI,:),:),1),[1,1]*length(CL))/size(T,2); 471CC.MDA.ACC00 = sum(CC.MDA.mmx(:,1:length(CL)+1:end),2)/sum(~isnan(cl)); 472CC.MDA.KAP00 = zeros(size(MDIX,1),1); 473 474CC.MD3.CMX00 = reshape(sum(CC.MD3.mmx(T(CC.TI,:),:),1),[1,1]*length(CL))/size(T,2); 475CC.MD3.ACC00 = sum(CC.MD3.mmx(:,1:length(CL)+1:end),2)/sum(~isnan(cl)); 476CC.MD3.KAP00 = zeros(size(GDIX,1),1); 477 478CC.LD2.CMX00 = reshape(sum(CC.LD2.mmx(T(CC.TI,:),:),1),[1,1]*length(CL))/size(T,2); 479CC.LD2.ACC00 = sum(CC.LD2.mmx(:,1:length(CL)+1:end),2)/sum(~isnan(cl)); 480CC.LD2.KAP00 = zeros(size(LD2IX,1),1); 481 482CC.LD3.CMX00 = reshape(sum(CC.LD3.mmx(T(CC.TI,:),:),1),[1,1]*length(CL))/size(T,2); 483CC.LD3.ACC00 = sum(CC.LD3.mmx(:,1:length(CL)+1:end),2)/sum(~isnan(cl)); 484CC.LD3.KAP00 = zeros(size(LD3IX,1),1); 485 486CC.LD4.CMX00 = reshape(sum(CC.LD4.mmx(T(CC.TI,:),:),1),[1,1]*length(CL))/size(T,2); 487CC.LD4.ACC00 = sum(CC.LD4.mmx(:,1:length(CL)+1:end),2)/sum(~isnan(cl)); 488CC.LD4.KAP00 = zeros(size(LD4IX,1),1); 489 490CC.LLH.CMX00 = reshape(sum(CC.LLH.mmx(T(CC.TI,:),:),1),[1,1]*length(CL))/size(T,2); 491CC.LLH.ACC00 = sum(CC.LLH.mmx(:,1:length(CL)+1:end),2)/sum(~isnan(cl)); 492CC.LLH.KAP00 = zeros(size(LLIX,1),1); 493for k = 1:size(MDIX,1), 494 CC.MDA.KAP00(k) = kappa(reshape(CC.MDA.mmx(k,:),[1,1]*length(CL))); 495 CC.MD3.KAP00(k) = kappa(reshape(CC.MD3.mmx(k,:),[1,1]*length(CL))); 496 CC.LD2.KAP00(k) = kappa(reshape(CC.LD2.mmx(k,:),[1,1]*length(CL))); 497 CC.LD3.KAP00(k) = kappa(reshape(CC.LD3.mmx(k,:),[1,1]*length(CL))); 498 CC.LD4.KAP00(k) = kappa(reshape(CC.LD4.mmx(k,:),[1,1]*length(CL))); 499 CC.LLH.KAP00(k) = kappa(reshape(CC.LLH.mmx(k,:),[1,1]*length(CL))); 500end; 501 502if length(CL) > 2, 503 return; 504end; 505 506 507if bitand(SWITCH,1), 508 CC.LLH.ERR00 = (mean(sign(JKLL),2)+1)/2; 509 CC.LDA.ERR00 = (mean(sign(JKLD),2)+1)/2; 510 CC.MD3.ERR00 = (mean(sign(JKGD),2)+1)/2; 511 CC.MDA.ERR00 = (mean(sign(JKD1-JKD2),2)+1)/2; 512 CC.GRB.ERR00 = (mean(sign(exp(-JKD2/2)-exp(-JKD1/2)),2)+1)/2; 513end; 514 515 516d=JKLD; 517tmp1 = d(1-min(T(:))+T(CC.TI,:),cl==CL(1)); 518[sum0,n0,ssq0] = sumskipnan(tmp1(:)); 519tmp2 = d(1-min(T(:))+T(CC.TI,:),cl==CL(2)); 520[sum1,n1,ssq1] = sumskipnan(tmp2(:)); 521CC.LDA.AUC = auc(tmp1,tmp2); 522CC.LDA.ERR(1,1) = mean(sign([tmp1(:)]))/2+1/2; 523CC.LDA.ERR(1,2) = mean(sign([tmp2(:)]))/2+1/2; 524CC.LDA.ERR(2,1) = mean(sign([mean(tmp1,1)']))/2+1/2; 525CC.LDA.ERR(2,2) = mean(sign([mean(tmp2,1)']))/2+1/2; 526s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 527s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 528s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 529SNR = 2*s./(s0+s1); % this is SNR+1 530CC.LDA.I = log2(SNR)/2; 531CC.LDA.SNR = SNR - 1; 532if 0, 533 clear tmp1 tmp2; 534 tmp1 = stat2(d(:,cl==CL(1)),2); 535 tmp2 = stat2(d(:,cl==CL(2)),2); 536 CC.LDA.TSD=stat2res(tmp1,tmp2); 537 CC.LDA.TSD.ERR=1/2-mean(sign([-d(:,cl==CL(1)),d(:,cl==CL(2))]),2)/2; 538elseif bitand(SWITCH,1), 539 CC.LDA.TSD=bci3eval(d(:,cl==CL(1)),d(:,cl==CL(2)),2); 540end; 541 542d = JKD1 - JKD2; 543tmp1 = d(1-min(T(:))+T(CC.TI,:),cl==CL(1)); 544[sum0,n0,ssq0] = sumskipnan(tmp1(:)); 545tmp2 = d(1-min(T(:))+T(CC.TI,:),cl==CL(2)); 546[sum1,n1,ssq1] = sumskipnan(tmp2(:)); 547CC.MDA.AUC = auc(tmp1,tmp2); 548CC.MDA.ERR(1,1) = mean(sign([tmp1(:)]))/2+1/2; 549CC.MDA.ERR(1,2) = mean(sign([tmp2(:)]))/2+1/2; 550CC.MDA.ERR(2,1) = mean(sign([mean(tmp1,1)']))/2+1/2; 551CC.MDA.ERR(2,2) = mean(sign([mean(tmp2,1)']))/2+1/2; 552s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 553s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 554s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 555SNR = 2*s./(s0+s1); % this is SNR+1 556CC.MDA.I = log2(SNR)/2; 557CC.MDA.SNR = SNR - 1; 558if 0, 559 clear tmp1 tmp2; 560 tmp1 = stat2(d(:,cl==CL(1)),2); 561 tmp2 = stat2(d(:,cl==CL(2)),2); 562 CC.MDA.TSD=stat2res(tmp1,tmp2); 563 CC.MDA.TSD.ERR=1/2-mean(sign([-d(:,cl==CL(1)),d(:,cl==CL(2))]),2)/2; 564elseif bitand(SWITCH,1), 565 CC.MDA.TSD=bci3eval(d(:,cl==CL(1)),d(:,cl==CL(2)),2); 566end; 567 568d = JKD6 - JKD5; 569tmp1 = d(1-min(T(:))+T(CC.TI,:),cl==CL(1)); 570[sum0,n0,ssq0] = sumskipnan(tmp1(:)); 571tmp2 = d(1-min(T(:))+T(CC.TI,:),cl==CL(2)); 572[sum1,n1,ssq1] = sumskipnan(tmp2(:)); 573CC.MD3.AUC = auc(tmp1,tmp2); 574CC.MD3.ERR(1,1) = mean(sign([tmp1(:)]))/2+1/2; 575CC.MD3.ERR(1,2) = mean(sign([tmp2(:)]))/2+1/2; 576CC.MD3.ERR(2,1) = mean(sign([mean(tmp1,1)']))/2+1/2; 577CC.MD3.ERR(2,2) = mean(sign([mean(tmp2,1)']))/2+1/2; 578s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 579s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 580s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 581SNR = 2*s./(s0+s1); % this is SNR+1 582CC.MD3.I = log2(SNR)/2; 583CC.MD3.SNR = SNR - 1; 584if 0, 585 clear tmp1 tmp2; 586 tmp1 = stat2(d(:,cl==CL(1)),2); 587 tmp2 = stat2(d(:,cl==CL(2)),2); 588 CC.MD3.TSD=stat2res(tmp1,tmp2); 589 CC.MD3.TSD.ERR=1/2-mean(sign([-d(:,cl==CL(1)),d(:,cl==CL(2))]),2)/2; 590elseif bitand(SWITCH,1), 591 CC.MD3.TSD=bci3eval(d(:,cl==CL(1)),d(:,cl==CL(2)),2); 592end; 593 594d = sqrt(JKD1) - sqrt(JKD2); 595tmp1 = d(1-min(T(:))+T(CC.TI,:),cl==CL(1)); 596[sum0,n0,ssq0] = sumskipnan(tmp1(:)); 597tmp2 = d(1-min(T(:))+T(CC.TI,:),cl==CL(2)); 598[sum1,n1,ssq1] = sumskipnan(tmp2(:)); 599CC.MD2.AUC = auc(tmp1,tmp2); 600CC.MD2.ERR(1,1) = mean(sign([tmp1(:)]))/2+1/2; 601CC.MD2.ERR(1,2) = mean(sign([tmp2(:)]))/2+1/2; 602CC.MD2.ERR(2,1) = mean(sign([mean(tmp1,1)']))/2+1/2; 603CC.MD2.ERR(2,2) = mean(sign([mean(tmp2,1)']))/2+1/2; 604s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 605s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 606s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 607SNR = 2*s./(s0+s1); % this is SNR+1 608CC.MD2.I = log2(SNR)/2; 609CC.MD2.SNR = SNR - 1; 610if 0, 611 clear tmp1 tmp2; 612 tmp1 = stat2(d(:,cl==CL(1)),2); 613 tmp2 = stat2(d(:,cl==CL(2)),2); 614 CC.MD2.TSD=stat2res(tmp1,tmp2); 615 CC.MD2.TSD.ERR=1/2-mean(sign([-d(:,cl==CL(1)),d(:,cl==CL(2))]),2)/2; 616elseif bitand(SWITCH,1), 617 CC.MD2.TSD=bci3eval(d(:,cl==CL(1)),d(:,cl==CL(2)),2); 618end; 619 620%%% 621if any(isnan(cl)), 622 t = perm(TRIG(isnan(cl)),T(CC.TI,:)); 623 t = t(t<=size(D,1)); 624 tmp= rs(D(t(:),:),size(T,2),1); 625 [CC.OUT.LDA] = ldbc(CC.MD,tmp); 626 CC.OUT.LDAcl = CL((CC.OUT.LDA>0)+1); 627 [CC.OUT.MDA] = mdbc(CC.MD,tmp); 628 [tmp,ix] = min(CC.OUT.MDA,[],2); 629 tmp = isnan(tmp); 630 ix(tmp) = NaN; % invalid output, not classified 631 ix(~tmp) = CL(ix(~tmp)); 632 CC.OUT.MDAcl = ix; 633end; 634 635return; 636 637d = JKD3 - JKD4; 638tmp1 = d(1-min(T(:))+T(CC.TI,:),cl==CL(1)); 639[sum0,n0,ssq0] = sumskipnan(tmp1(:)); 640tmp2 = d(1-min(T(:))+T(CC.TI,:),cl==CL(2)); 641[sum1,n1,ssq1] = sumskipnan(tmp2(:)); 642CC.MLL.AUC = auc(tmp1,tmp2); 643CC.MLL.ERR(1,1) = mean(sign([tmp1(:)]))/2+1/2; 644CC.MLL.ERR(1,2) = mean(sign([tmp2(:)]))/2+1/2; 645CC.MLL.ERR(2,1) = mean(sign([mean(tmp1,1)']))/2+1/2; 646CC.MLL.ERR(2,2) = mean(sign([mean(tmp2,1)']))/2+1/2; 647s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 648s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 649s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 650SNR = 2*s./(s0+s1); % this is SNR+1 651CC.MLL.I = log2(SNR)/2; 652CC.MLL.SNR = SNR - 1; 653if 0, 654 clear tmp1 tmp2; 655 tmp1 = stat2(d(:,cl==CL(1)),2); 656 tmp2 = stat2(d(:,cl==CL(2)),2); 657 CC.MLL.TSD=stat2res(tmp1,tmp2); 658 CC.MLL.TSD.ERR=mean(sign([-d(:,cl==CL(1)),d(:,cl==CL(2))]),2)/2+1/2; 659elseif bitand(SWITCH,1), 660 CC.MLL.TSD=bci3eval(d(:,cl==CL(1)),d(:,cl==CL(2)),2); 661end; 662 663d = exp(-JKD1/2)-exp(-JKD2/2); 664tmp1 = d(1-min(T(:))+T(CC.TI,:),cl==CL(1)); 665[sum0,n0,ssq0] = sumskipnan(tmp1(:)); 666tmp2 = d(1-min(T(:))+T(CC.TI,:),cl==CL(2)); 667[sum1,n1,ssq1] = sumskipnan(tmp2(:)); 668CC.GRB.AUC = auc(tmp1,tmp2); 669CC.GRB.ERR(1,1) = mean(sign([tmp1(:)]))/2+1/2; 670CC.GRB.ERR(1,2) = mean(sign([tmp2(:)]))/2+1/2; 671CC.GRB.ERR(2,1) = mean(sign([mean(tmp1,1)']))/2+1/2; 672CC.GRB.ERR(2,2) = mean(sign([mean(tmp2,1)']))/2+1/2; 673s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 674s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 675s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 676SNR = 2*s./(s0+s1); % this is SNR+1 677CC.GRB.I = log2(SNR)/2; 678CC.GRB.SNR = SNR - 1; 679if 0, 680 clear tmp1 tmp2; 681 tmp1 = stat2(d(:,cl==CL(1)),2); 682 tmp2 = stat2(d(:,cl==CL(2)),2); 683 CC.GRB.TSD=stat2res(tmp1,tmp2); 684 CC.GRB.TSD.ERR=1/2-mean(sign([-d(:,cl==CL(1)),d(:,cl==CL(2))]),2)/2; 685elseif bitand(SWITCH,1), 686 CC.GRB.TSD=bci3eval(d(:,cl==CL(1)),d(:,cl==CL(2)),2); 687end; 688 689