1function [CC,Q,tsd,md,cc]=fc0(D,TRIG,T,arg4) 2% FC finds a classifier for asnychroneous data 3% 4% [CC,Q,TSD,MD]=fc0(D,TRIG,class_times [,SWITCH]); 5% 6% D data, each row is one time point 7% TRIG trigger time points 8% class_times each row defines a segment used for Classification 9% SWITCH 0 [default] minimizes result 10% 2 only have of the possible t1-t2 combinations are used 11% 12% CC contains LDA and MD classifiers 13% Q is a list of classification quality for each time of 'class_times' 14% TSD returns the LDA classification 15% MD returns the MD classification 16% 17% 18% [CC,Q,TSD,MD]=fc(AR,find(trig>0.5)-257,reshape(1:14*128,16,14*8)'); 19% 20% see also: COVM.M, QCMAHAL, MDBC, LDBC, 21 22% 23% $Revision: 1.1 $ 24% $Id$ 25% Copyright (c) 1999-2002, 2004 by Alois Schloegl <alois.schloegl@gmail.com> 26% This is part of the BIOSIG-toolbox http://biosig.sf.net/ 27% 28 29% CHANGELOG 30% 10.12.2001 changed Covariance-Matrix from unbiased to biased 31% 20.12.2001 ROC and AUC included 32% 04.01.2002 zeros class added 33% md changed to diff of log of MD 34% 18.01.2002 Gaussian Radial basis functions 35% 12.02.2002 included in BCI7-paradigm 36% 13.02.2002 CC.D included 37% 30.04.2002 JACKKNIFE included 38% 19.07.2002 SWITCH included 39 40% This program is free software; you can redistribute it and/or 41% modify it under the terms of the GNU General Public License 42% as published by the Free Software Foundation; either version 2 43% of the License, or (at your option) any later version. 44% 45% This program is distributed in the hope that it will be useful, 46% but WITHOUT ANY WARRANTY; without even the implied warranty of 47% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 48% GNU General Public License for more details. 49% 50% You should have received a copy of the GNU General Public License 51% along with this program; if not, write to the Free Software 52% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 53 54if nargin>3; 55 SWITCH = arg4; 56end; 57 58TRIG=TRIG(:); 59if ~all(D(:,1)==1) 60% D1=[ones(size(D,1)-1,1),diff(D)]; 61 D =[ones(size(D,1),1),D]; 62%else 63% D1=[ones(size(D,1)-1,1),diff(D(:,2:end))]; 64end; 65 66dT = (min(T(:)):max(T(:))); 67nc = length(dT); 68%%% pad sufficient NaN's 69off=min(TRIG(1)+min(T(:))-1,0); 70D=[repmat(nan,-off,size(D,2));D]; 71TRIG=TRIG-off; 72off=max(max(T(:))+TRIG(length(TRIG))-length(D),0); 73D=[D;repmat(nan,off,size(D,2))]; 74 75 76[CC0,NN0] = covm(D,'E'); 77for k = 1:size(T,1), 78 t = perm(TRIG,T(k,:)); 79 t = t(:); 80 [cc{k},nn{k}] = covm(D(t,:),'M'); 81 %CC{k}=cc{k}./nn{k}; 82end; 83[Q , d] = qcmahal(cc); 84if bitand(SWITCH,1) 85 CC.d=d; 86end; 87if bitand(SWITCH,2), 88 d=d+d'; 89 d=d(:,1:size(d,2)/2); 90end; 91 92[ix,iy] = find(d == max(d(:))); 93 94ix = ix(1); iy = iy(1); 95if ix>iy, 96 tmp=ix;ix=iy;iy=tmp; 97end; 98CC.TI = [ix,iy]; 99CC.MD = {cc{ix},cc{iy}}; 100[CC.Q,CC.D]= qcmahal(CC.MD); 101Q=CC.Q; 102 103 104%CC.IR = mdbc(CC.MD); 105 106JKD1=zeros(length(dT),length(TRIG)); 107JKD2=zeros(length(dT),length(TRIG)); 108%CC.SGN = zeros(length(TRIG),nc); 109%% jackknife 110for l=1:length(TRIG), 111 T1 = TRIG(l) + T(CC.TI(1),:); 112 T2 = TRIG(l) + T(CC.TI(2),:); 113 114 tmp = D(T1,:); 115 [tmp1,tmp2]=covm(tmp,'M'); 116 cc1 = CC.MD{1}-tmp1; 117 118 tmp = D(T2,:); 119 [tmp1,tmp2]=covm(tmp,'M'); 120 cc2 = CC.MD{2}-tmp1; 121 122 t = TRIG(l) + dT; 123 d = mdbc({cc1,cc2},D(t,:)); 124 JKD1(:,l)=d(:,1); 125 JKD2(:,l)=d(:,2); 126 127 d = llbc({cc1,cc2},D(t,:)); 128 JKD3(:,l)=d(:,1); 129 JKD4(:,l)=d(:,2); 130 131 JKLD(:,l) = ldbc({cc1,cc2}, D(t,:)); 132end; 133if bitand(SWITCH,1), 134 CC.LDA.ERR00 = (mean(sign(JKLD),2)+1)/2; 135 CC.MDA.ERR00 = (mean(sign(JKD1-JKD2),2)+1)/2; 136 CC.GRB.ERR00 = (mean(sign(exp(-JKD2/2)-exp(-JKD1/2)),2)+1)/2; 137end; 138 139d = JKD3 - JKD4; 140tmp1 = d(1-min(T(:))+T(CC.TI(1),:),:); 141[sum0,n0,ssq0] = sumskipnan(tmp1(:)); 142tmp2 = d(1-min(T(:))+T(CC.TI(2),:),:); 143[sum1,n1,ssq1] = sumskipnan(tmp2(:)); 144CC.MLL.AUC = auc(tmp1,tmp2); 145CC.MLL.ERR(1,:) = mean(sign([tmp1(:),tmp2(:)]))/2+1/2; 146CC.MLL.ERR(2,:) = mean(sign([mean(tmp1,1)',mean(tmp2,1)']))/2+1/2; 147s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 148s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 149s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 150SNR = 2*s./(s0+s1); % this is SNR+1 151CC.MLL.I = log2(SNR)/2; 152CC.MLL.SNR = SNR - 1; 153clear tmp; tmp.datatype='STAT2'; 154[tmp.SUM,tmp.N,tmp.SSQ] = sumskipnan(d,2); 155if bitand(SWITCH,1), 156 CC.MLL.TSD=tmp; 157end; 158 159d = JKD1 - JKD2; 160tmp1 = d(1-min(T(:))+T(CC.TI(1),:),:); 161[sum0,n0,ssq0] = sumskipnan(tmp1(:)); 162tmp2 = d(1-min(T(:))+T(CC.TI(2),:),:); 163[sum1,n1,ssq1] = sumskipnan(tmp2(:)); 164CC.MDA.AUC = auc(tmp1,tmp2); 165CC.MDA.ERR(1,:) = mean(sign([tmp1(:),tmp2(:)]))/2+1/2; 166CC.MDA.ERR(2,:) = mean(sign([mean(tmp1,1)',mean(tmp2,1)']))/2+1/2; 167s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 168s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 169s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 170SNR = 2*s./(s0+s1); % this is SNR+1 171CC.MDA.I = log2(SNR)/2; 172CC.MDA.SNR = SNR - 1; 173clear tmp; tmp.datatype='STAT2'; 174[tmp.SUM,tmp.N,tmp.SSQ] = sumskipnan(d,2); 175if bitand(SWITCH,1), 176 CC.MDA.TSD=tmp; 177end; 178 179d = exp(-JKD1/2)-exp(-JKD2/2); 180tmp1 = d(1-min(T(:))+T(CC.TI(1),:),:); 181[sum0,n0,ssq0] = sumskipnan(tmp1(:)); 182tmp2 = d(1-min(T(:))+T(CC.TI(2),:),:); 183[sum1,n1,ssq1] = sumskipnan(tmp2(:)); 184CC.GRB.AUC = auc(tmp1,tmp2); 185CC.GRB.ERR(1,:) = mean(sign([tmp1(:),tmp2(:)]))/2+1/2; 186CC.GRB.ERR(2,:) = mean(sign([mean(tmp1,1)',mean(tmp2,1)']))/2+1/2; 187s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 188s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 189s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 190SNR = 2*s./(s0+s1); % this is SNR+1 191CC.GRB.I = log2(SNR)/2; 192CC.GRB.SNR = SNR - 1; 193clear tmp; tmp.datatype='STAT2'; 194[tmp.SUM,tmp.N,tmp.SSQ] = sumskipnan(d,2); 195if bitand(SWITCH,1), 196 CC.GRB.TSD=tmp; 197end; 198 199d=JKLD; 200tmp1 = d(1-min(T(:))+T(CC.TI(1),:),:); 201[sum0,n0,ssq0] = sumskipnan(tmp1(:)); 202tmp2 = d(1-min(T(:))+T(CC.TI(2),:),:); 203[sum1,n1,ssq1] = sumskipnan(tmp2(:)); 204CC.LDA.AUC = auc(tmp1,tmp2); 205CC.LDA.ERR(1,:) = mean(sign([tmp1(:),tmp2(:)]))/2+1/2; 206CC.LDA.ERR(2,:) = mean(sign([mean(tmp1,1)',mean(tmp2,1)']))/2+1/2; 207s0 = (ssq0-sum0.*sum0./n0)./(n0-1); 208s1 = (ssq1-sum1.*sum1./n1)./(n1-1); 209s = (ssq0+ssq1-(sum0+sum1).*(sum0+sum1)./(n0+n1))./(n0+n1-1); 210SNR = 2*s./(s0+s1); % this is SNR+1 211CC.LDA.I = log2(SNR)/2; 212CC.LDA.SNR = SNR - 1; 213clear tmp; tmp.datatype='STAT2'; 214[tmp.SUM,tmp.N,tmp.SSQ] = sumskipnan(d,2); 215if bitand(SWITCH,1), 216 CC.LDA.TSD=tmp; 217end; 218 219 220 221 222return 223 224%%%% reference intervall 225if nargin>4, 226 T = perm(TRIG,t(k,:)); 227 T = T(T<=size(D,1)); 228 tmp = D(T(:),:); 229 C0r = tmp'*tmp; 230end; 231 232for k = 1:size(t,1), 233 for l = 1:length(CL), 234 T = perm(TRIG(cl==CL(l)),t(k,:)); 235 T = T(T<=size(D,1)); 236 tmp = D(T(:),:); 237 C{k,l} = tmp'*tmp; 238 end; 239 %[Q(k),d{k}] = qcmahal({C0r,C{k,:}}); 240 [Q(k),d{k}] = qcmahal({C{k,:}}); 241 lnQ(k) = mean(log(d{k}(~eye(length(d{k}))))); 242end; 243[maxQ,CC.TI] = max(Q); %d{K}, 244%CC.TI = K; 245CC.MD = {C{CC.TI,:}}; 246CC.IR = mdbc({C{CC.TI,:}}); 247CC.D = d{CC.TI}; 248Q = Q(CC.TI); 249 250[maxQ,CC.lnTI] = max(lnQ); %d{K}, 251CC.DistMXln = d{CC.lnTI}; 252CC.MDln = {C{CC.lnTI,:}}; 253 254% LDA 255C0 = zeros(size(C{CC.TI,1})); 256for l=1:length(CL); 257 [M{l},sd,COV,xc,N,R2] = decovm(C{CC.TI,l}); 258 C0 = C0 + C{CC.TI,l}; 259end; 260[M0,sd,COV0,xc,N,R2] = decovm(C0); 261w = COV0\(M{1}'-M{2}'); 262w0 = M0*w; 263CC.LDA.b = w0; 264CC.LDA.w = -w; 265CC.lda = [w0; -w]; 266 267% MD 268md = mdbc(CC.MD,D); 269 270 271[tmp,IX] = min(md,[],2); 272for k = 1:size(t,1), 273 H0{k} = zeros(length(CL)); 274 for l = 1:size(t,2), 275 T = TRIG+t(k,l); 276 T = T(T<=size(D,1)); 277 [tmp,sd,H] = kappa(cl(:)-min(cl)+1,IX(T),length(CL)); 278 H0{k} = H0{k} + H; 279 end; 280 [CC.KAPPA(k),se,tmp,z,CC.OA(k),CC.SA(k,1:2)] = kappa(H0{k}); 281end; 282CC.H0=H0; 283 284if length(CL)>2, return; end; 285 286CC.LDA=eval_offline(D*CC.lda,TRIG,cl); 287CC.MDA=eval_offline(md(:,1)-md(:,2),TRIG,cl); 288CC.GRB=eval_offline(exp(-md(:,2)/2)-exp(-md(:,1)/2),TRIG,cl); 289CC.GR1=eval_offline(exp(-md(:,2))-exp(-md(:,1)),TRIG,cl); 290 291nc = ceil(max(diff(TRIG))*1.5); 292ERR00md = zeros(1,nc); 293ERR00gr = zeros(1,nc); 294ERR00ld = zeros(1,nc); 295N00 = zeros(1,nc); 296 297%% jackknife 298for l=1:length(cl), 299 c = find(cl(l)==CL); 300 T = TRIG(l)+t(CC.TI,:); 301 T = T(T<=size(D,1)); 302 tmp = D(T(:),:); 303 cc{c} = CC.MD{c}-tmp'*tmp; 304 cc{3-c} = CC.MD{3-c}; 305 T = TRIG(l)+(1:nc); 306 T = T(T<=size(D,1)); 307 [tmp] = mdbc(cc,D(T,:)); 308 ERR00md(1:length(T))= ERR00md(1:length(T))+sign(1.5-c)*sign(tmp(:,2)-tmp(:,1))'; 309 ERR00gr(1:length(T))= ERR00gr(1:length(T))+sign(1.5-c)*sign(exp(-tmp(:,1)/2)-exp(-tmp(:,2)/2))'; 310 tmp = ldbc(cc,D(T,:)); 311 ERR00ld(1:length(T))= ERR00ld(1:length(T))+sign(c-1.5)*sign(tmp)'; 312 N00(1:length(T)) = N00(1:length(T))+1; 313end; 314 315CC.ERR00md = (1-ERR00md./N00)/2; 316CC.ERR00gr = (1-ERR00gr./N00)/2; 317CC.ERR00ld = (1-ERR00ld./N00)/2; 318CC.N00 = N00; 319 320