findclassifier2.m

matlab数字信号处理工具箱

function [CC,Q,tsd,md]=findclassifier2(D,TRIG,cl,T,t0,SWITCH)
% FINDCLASSIFIER2
%       is very similar to FINDCLASSIFIER1 but uses a different 
%       criterion for selecting the time segment.
% [CC,Q,TSD,MD]=findclassifier2(D,TRIG,Class,class_times,t_ref);
%
% D 	data, each row is one time point
% TRIG	trigger time points
% Class class information
% class_times	classification times, combinations of times must be in one row 
% t_ref	reference time for Class 0 (optional)
%
% CC 	contains LDA and MD classifiers
% Q  	is a list of classification quality for each time of 'class_times'
% TSD 	returns the LDA classification 
% MD	returns the MD  classification 
%
% [CC,Q,TSD,MD]=findclassifier2(AR,find(trig>0.5)-257,~mod(1:80,2),reshape(1:14*128,16,14*8)');
%
%
% Reference(s): 
% [1] Schl鰃l A., Neuper C. Pfurtscheller G.
% 	Estimating the mutual information of an EEG-based Brain-Computer-Interface
%  	Biomedizinische Technik 47(1-2): 3-8, 2002.
% [2] A. Schl鰃l, C. Keinrath, R. Scherer, G. Pfurtscheller,
%	Information transfer of an EEG-based Bran-computer interface.
%	Proceedings of the 1st International IEEE EMBS Conference on Neural Engineering, Capri, Italy, Mar 20-22, 2003 


%   Copyright (C) 1999-2004 by Alois Schloegl <[e-mail]>	
%	$Id: findclassifier2.m,v 1.3 2004/09/02 22:12:14 schloegl Exp $


% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the  License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

tsd=[];md=[];

if nargin<6,
        SWITCH=0;
end;
if nargin>4,
        if isempty(t0),
                t0=logical(ones(size(T,1),1));
        end;
end;
tmp=cl;tmp(isnan(tmp))=0;
if any(rem(tmp,1) & ~isnan(cl)),
        fprintf(2,'Error %s: class information is not integer\n',mfilename);
        return;
end;
if length(TRIG)~=length(cl);
        fprintf(2,'number of Triggers do not match class information');
end;

CL = unique(cl(~isnan(cl)));

CL = sort(CL);
TRIG = TRIG(:);
if ~all(D(:,1)==1)
        %        D1=[ones(size(D,1)-1,1),diff(D)];
        D = [ones(size(D,1),1),D];
        %else
        %	D1=[ones(size(D,1)-1,1),diff(D(:,2:end))];
end;

% add sufficient NaNs at the beginning and the end 
tmp = min(TRIG)+min(min(T))-1;
if tmp<0,
        TRIG = TRIG - tmp;
        D = [repmat(nan,[-tmp,size(D,2)]);D];
end;        
tmp = max(TRIG)+max(max(T))-size(D,1);
if tmp>0,
        D = [D;repmat(nan,[tmp,size(D,2)])];
end;        

% estimate classification result for all time segments in T - without crossvalidation 
CMX = zeros([size(T,1),length(CL)*[1,1]]);
for k = 1:size(T,1),
        cmx = zeros(length(CL));
        for l = 1:length(CL), 
                t = perm(TRIG(cl==CL(l)),T(k,:));
                %t = t(t<=size(D,1));
                [C{k,l},tmp] = covm(D(t(:),:),'M');
        end;
        %[Q(k),d{k}] = qcmahal({C0r,C{k,:}});
        [CC.QC(k),d{k}] = qcmahal({C{k,:}});
        %[CC.QC3(k,1:2),d3{k}] = qcloglik({C{k,:}});
        CC.QC2(k) = sum(sqrt(d{k}(:)))/(size(d{k},1)*(size(d{k},1)-1));
        lnQ(k) = mean(log(d{k}(~eye(length(d{k})))));
        for l = 1:length(CL), 
                t = perm(TRIG(cl==CL(l)),T(k,:));
                %t = t(t<=size(D,1));
                [tmp] = mdbc({C{k,:}},D(t(:),:));
                [tmp,ix] = min(tmp,[],2);
                tmp = isnan(tmp);
                ix(tmp) = NaN; %NC(1)+1;	% not classified; any value but not 1:length(MD)
                ix(~tmp) = CL(ix(~tmp));
                tmp = histo3([ix;CL(:)]);
                cmx(tmp.X,l) = tmp.H-1;            
        end;
        CMX(k,:,:)  = cmx;
        CC.KAPPA(k) = kappa(cmx);
        CC.ACC(k)   = sum(diag(cmx))/sum(cmx(:));

        t = perm(TRIG,T(k,:));
        tmp = repmat(cl(:)',size(T,2),1);
        [tmp,acc_g(k,:),kap_g(k,:),H] = gdbc({C{k,:}},D(t(:),:),tmp(:));
end;	
CC.GDBC.acc = acc_g;
CC.GDBC.kap = kap_g;
CC.GDBC.H = H;
CC.GDBC.p = tmp;

% identify best classification time 
if nargin>4,
        tmp = CC.QC;
        tmp(~t0) = 0;
        [maxQ,CC.TI] = max(tmp); %d{K},
else
        [maxQ,CC.TI] = max(CC.QC); %d{K},
end;
CC.qc = [CC.QC',CC.KAPPA',CC.ACC',CC.QC2'];
%plot([CC.QC',CC.KAPPA',CC.ACC',CC.QC2',CC.QC3',t0]);drawnow
[maxQ(1),TI(1)] = max(CC.QC'.*t0); %d{K},
[maxQ(2),TI(2)] = max(CC.KAPPA'.*t0); %d{K},
[maxQ(3),TI(3)] = max(CC.ACC'.*t0); %d{K},
[maxQ(4),TI(4)] = max(CC.QC2'.*t0); %d{K},
%[maxQ(5),TI(5)] = max(CC.QC3(:,1).*t0); %d{K},
%[maxQ(6),TI(6)] = max(CC.QC3(:,2).*t0); %d{K},
CC.TIS = TI;
CC.TI  = TI(2);
%if any(TI(1)~=TI),
%        fprintf(1,'FINDCLASSIFIER2: \nTIS:  %i\t%i\t%i\t%i',TI);
%        fprintf(1,'\n     %5.3f\t%5.3f\t%5.3f\t%5.3f\n',maxQ);
%end;

% build MD classifier 
CC.MD  = {C{CC.TI,:}};
CC.IR  = mdbc({C{CC.TI,:}});
CC.D   = d{CC.TI};
CC.Q   = CC.QC(CC.TI);
CC.CMX = squeeze(CMX(CC.TI,:,:));
CC.CMX = CMX;
CC.tsc = T(CC.TI,[1,end]);

m1=decovm(CC.MD{1});
m2=decovm(CC.MD{2});
tmp=mdbc(CC.MD,[1,m1;1,m2]);
CC.scale=[1,1]*max(abs(tmp(:)));  % element 1 

[maxQ,CC.lnTI] = max(lnQ); %d{K},
CC.DistMXln = d{CC.lnTI};
CC.MDln = {C{CC.lnTI,:}};

% alternative classifier using two different time segments.
if 1,
        [Q,d] = qcmahal({C{:}}');
        CC.T2.D = d;
        [ix,iy] = find(d==max(d(:)));
        ix=mod(ix(1)-1,size(C,1))+1;
        iy=mod(iy(1)-1,size(C,1))+1;
        CC.T2.TI = [ix(1),iy(1)];
        CC.T2.MD = {C{ix(1),1},C{iy(1),2}};
end;

% build LDA classifier
if length(CL)==2,
        % LDA 
        C0 = zeros(size(C{CC.TI,1}));
        for l=1:length(CL);
                [M{l},sd,COV,xc,N,R2] = decovm(C{CC.TI,l});
                C0 = C0 + C{CC.TI,l};
        end;
        [M0,sd,COV0,xc,N,R2] = decovm(C0);
        w     = COV0\(M{1}'-M{2}');
        w0    = M0*w;
        %CC.LDA.b = w0;
        %CC.LDA.w = -w;
        CC.lda = [w0; -w];
        
        % MD 
        tsd = ldbc(CC.MD,D);
end;
md = mdbc(CC.MD,D);
lld= llbc(CC.MD,D);

% bias correction not used anymore.
CC.BIAS.LDA=0;%mean(tsd);
CC.BIAS.MDA=0;%mean(diff(-md,[],2));
CC.BIAS.GRB=0;%mean(diff(-exp(-md/2),[],2));

% define datatype
CC.datatype = 'Classifier';
CC.Classes = CL;


%% cross-validation with jackknife (leave one trial out)
nc  = max(max(T))-min(min(T))+1;

JKD   = repmat(nan,[nc,length(CL),length(TRIG)]);
JKLL  = JKD;
JKD1  = repmat(nan,[nc,length(TRIG)]);
JKD2  = repmat(nan,[nc,length(TRIG)]);
JKD3  = repmat(nan,[nc,length(TRIG)]);
JKD4  = repmat(nan,[nc,length(TRIG)]);
JKD5  = repmat(nan,[nc,length(TRIG)]);
JKD6  = repmat(nan,[nc,length(TRIG)]);
JKLD  = repmat(nan,[nc,length(TRIG)]);
for l = find(~isnan(cl(:)'));1:length(cl);
        c = find(cl(l)==CL);
        t = TRIG(l)+T(CC.TI,:);
        %t = t(t<=size(D,1));
        [tmp,tmpn] = covm(D(t(:),:),'M');
        
        cc 	= CC.MD;
        cc{c}   = CC.MD{c}-tmp;
        
        %t = TRIG(l)+(1:nc);
        %t = t(t<=size(D,1));
        t = TRIG(l)+(min(min(T)):max(max(T)));
        
        d = llbc(cc,D(t,:));
        JKLL(:,:,l) = d;
        [tmp,LLIX(:,l)] = max(d,[],2);
        if length(CL)==2,
                JKD3(:,l)=d(:,1);
                JKD4(:,l)=d(:,2);
        end;
        
        d = mdbc(cc,D(t,:));
        JKD(:,:,l) = d;
        [tmp,MDIX(:,l)] = min(d,[],2);
        
        if length(CL)==2,
                JKD1(:,l) = d(:,1);
                JKD2(:,l) = d(:,2);
	end;
	                
        d = gdbc(cc,D(t,:));
        JKGD(:,:,l) = d;
        [tmp,GDIX(:,l)] = max(d,[],2);
        
        if length(CL)==2,
                JKD5(:,l) = d(:,1);
                JKD6(:,l) = d(:,2);
                
                LDA(:,l) = ldbc(cc);
                JKLD(:,l) = D(t,:)*LDA(:,l);
        end;	
end;
if length(CL)==2,
        [CC.ldaC0,NN] = covm(LDA','D0');
        %CC.ldaC0=CC.ldaC0./NN*min(0,sum(~isnan(CL))-1); 
        % since NN==min(0,sum(~isnan(CL))-1), no need to rescale
end;	

% Concordance matrix with cross-validation 
CC.lmx= zeros([size(LLIX,1),length(CL)^2]);
CC.LLH.I0 = zeros([size(LLIX,1),length(CL)]);
CC.LLH.I  = zeros([size(LLIX,1),1]);

CC.mmx= zeros([size(MDIX,1),length(CL)^2]);
CC.MD2.I0 = zeros([size(MDIX,1),length(CL)]);
CC.MD2.I  = zeros([size(MDIX,1),1]);
tmp   = zeros([size(MDIX,1),length(CL)]);
for k = 1:length(CL),

        jkd = squeeze(JKGD(:,k,:));
        o = bci3eval(jkd(:,cl~=k),jkd(:,cl==k),2);
        CC.MD3.TSD{k}  = o;
        CC.MD3.I0(:,k) = log2(2*var(jkd,[],2)./(var(jkd(:,cl==k),[],2) + var(jkd(:,cl~=k),[],2)))/2;
        [sum0,n0,ssq0]=sumskipnan(jkd(:,cl==k),2);
        [sum1,n1,ssq1]=sumskipnan(jkd(:,cl~=k),2);
        s0  = (ssq0-sum0.*sum0./n0)./(n0-1);
        s1  = (ssq1-sum1.*sum1./n1)./(n1-1);
        s   = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1);