% fisher linear discriminant    jp lewis 
% 
% Note this gives the line to project onto that maximally separates
% the data.  The separating hyperplane is perpendicular to this.
% The math:
% Maximize between-class scatter while holding total constant,
%
%    max(w)     w'(m1-m2)(m1-m2)'w  +  L(w'Sw - k)
%
% where L=lagrange multiplier,  k = constant.


hold off;
hold on;

% generate the data
N = 70;		% number of data points
data = zeros([N,3]);	% x,y, classification -1,1
for i=1:N
  x = 2*rand() - 1;
  y = 2*rand() - 1;
  data(i,1) = x;
  data(i,2) = y;
 
  % 'L' shaped classes
    if ((x < 0.) && (y < 0.8))
      data(i,3) = -1;
      plot(x,y,'o');
    else
      data(i,3) = 1;
      plot(x,y,'x');
    end
end

% means of data
  n1 = 0;		% number in class 1
  n2 = 0;		% class 2
  m1x = 0.; m1y = 0.;	% x,y mean of class 1
  m2x = 0.; m2y = 0.;	% x,y mean of class 2
  mx = 0;  my = 0;	% x,y mean of all

  for i=1:N
    mx = mx + data(i,1);
    my = my + data(i,2);

    if (data(i,3) == 1)
      m1x = m1x + data(i,1);
      m1y = m1y + data(i,2);
      n1 = n1 + 1;
    else
      m2x = m2x + data(i,1);
      m2y = m2y + data(i,2);
      n2 = n2 + 1;
    end
  end

  mx = mx / N;
  my = my / N;

  m1x = m1x / n1;
  m1y = m1y / n1;

  m2x = m2x / n2;
  m2y = m2y / n2;


% Between-class matrix
v = [m1x - m2x; m1y - m2y];
B = v*v';	% between-class variance

% total covariance matrix
T = zeros([2 2]);
u=[0;0];
for i=1:N
  u(1) = (data(i,1) - mx);
  u(2) = (data(i,2) - my);
  T = T + u*u';
end
T = (1./N) * T;
T

% Generalized Eigenvalue problem.  L=lagrange multiplier
% max(w):  w'Bw + L(w'Tw - 1)
% d/dw=0:  2Bw + 2LT2 = 0
% Bw = -LTw
% inv(T)Bw = -Lw
% result should be eigenvector corresponding to *largest* eigenvalue

M = inv(T)*B;		% convert to standard eigenvalue
[V,D] = eig(M);		
[W,D] = eig(B,T);	% generalized eigenvalue version

v = V(:,2);		% larger eigenvector is projection line,
w = W(:,1);		% smaller is orthogonal and is separating plane
plot([0,v(1)], [0, v(2)], 'r');		% plot the separating only
plot([0,w(1)], [0, w(2)], 'g');		% gen.eig solution is same