1 %solution_sheet12
   2 
   3 load features_lecture_adaptiveBCI
   4 
   5 %Ex 1
   6 icl1=find(fv_train.y(1,:));
   7 icl2=find(fv_train.y(2,:));
   8 icl1_test=find(fv_test.y(1,:));
   9 icl2_test=find(fv_test.y(2,:));
  10 
  11 %plot distributions of train and test datasets.
  12 plot(fv_train.x(1,icl1),fv_train.x(2,icl1),'rx')
  13 hold on
  14 plot(fv_train.x(1,icl2),fv_train.x(2,icl2),'go')
  15 plot(fv_test.x(1,icl1_test),fv_test.x(2,icl1_test),'kx')
  16 plot(fv_test.x(1,icl2_test),fv_test.x(2,icl2_test),'bo')
  17 
  18 %calculate LDA classifier of train data. new version
  19 mean_train(:,1)=mean(fv_train.x(:,icl1)')';
  20 mean_train(:,2)=mean(fv_train.x(:,icl2)')';
  21 
  22 w=pinv(cov(fv_train.x'));
  23 w=-w*(diff(mean_train,1,2));
  24 b=-mean(fv_train.x')*w;
  25 
  26 %calculate LDA classifier old version.
  27 w_old=pinv((cov(fv_train.x(:,icl1)')+cov(fv_train.x(:,icl2)'))/2)*(diff(mean_train,1,2));
  28 b_old=-mean(mean_train')*w_old;
  29 
  30 %both are equivalent
  31 w./w_old
  32 b/b_old
  33 
  34 %plot classifier boundary
  35 yl= (-w(1)*xlim-b)/w(2);
  36 plot(xlim, yl, 'k');
  37 yl= (-w_old(1)*xlim-b_old)/w_old(2);
  38 plot(xlim, yl, 'm');
  39 
  40 %calculate LDA error with train classifier.
  41 out=w'*fv_test.x+b*ones(1,size(fv_test.x,2));
  42 out_old=w_old'*fv_test.x+b_old*ones(1,size(fv_test.x,2));
  43 out_label=sign(out);
  44 err_orig=mean(out_label~=[-1 1]*fv_test.y);
  45 out_label=sign(out_old);
  46 err_old=mean(out_label~=[-1 1]*fv_test.y);
  47 
  48 %adapt trial by trial
  49 UC=0.01;
  50 %type one, using extended covariance matrix.
  51 [nf,N]=size(fv_train.x);
  52 fvold.x=[ones(1,N) ;fv_train.x];
  53 Cold=fvold.x*fvold.x'/N;
  54 invCold=pinv(Cold);
  55 invCnew=invCold;
  56 factor=1/(1-UC);
  57 out_test=zeros(1,size(fv_test.y,2));
  58 cadap.w=w;
  59 cadap.b=b;
  60  for itrial=1:size(fv_test.x,2)
  61      out_test(itrial)=cadap.w'*fv_test.x(:,itrial)+cadap.b;
  62      v=invCnew*[1 ;fv_test.x(:,itrial)];
  63      denominator=((1-UC)/UC)+[1;fv_test.x(:,itrial)]'*v;
  64      invCnew=factor*(invCnew-(v*v'/denominator));
  65      invCnew=(invCnew+invCnew')/2;
  66      cadap.w=invCnew(2:end,2:end)*diff(mean_train,1,2);
  67      cadap.b=invCnew(1,2:end)*diff(mean_train,1,2);
  68  end;
  69 
  70  yl= (-cadap.w(1)*xlim-cadap.b)/cadap.w(2);
  71 plot(xlim, yl, 'k');
  72 
  73 %calculate new errors
  74 out_label=sign(out_test);
  75 err= mean(out_label~=[-1 1]*fv_test.y);
  76 
  77  
  78  %second way: using normal covariance and updating the pooled mean
  79  
  80 icovnew=inv(cov(fv_train.x'));
  81 pmold=mean(fv_train.x');
  82 pmean=pmold';
  83 out_test=zeros(1,size(fv_test.y,2));
  84 cadap.w=w;
  85 cadap.b=b;
  86 for itrial=1:size(fv_test.x,2)
  87      out_test(itrial)=cadap.w'*fv_test.x(:,itrial)+cadap.b;
  88      pmean=(1-UC)*pmean+UC*fv_test.x(:,itrial);
  89      sample=fv_test.x(:,itrial)-pmean;
  90      v=icovnew*sample;
  91      denominator=((1-UC)/UC)+sample'*v;
  92      icovnew=factor*(icovnew-(v*v'/denominator));
  93      icovnew=(icovnew+icovnew')/2;
  94      cadap.w=icovnew*diff(mean_train,1,2);
  95      cadap.b=-cadap.w'*pmean;
  96  end
  97 
  98  yl= (-cadap.w(1)*xlim-cadap.b)/cadap.w(2);
  99 plot(xlim, yl, 'b');
 100 
 101 %calculte new errors
 102 out_label=sign(out_test);
 103 err= mean(out_label~=[-1 1]*fv_test.y);