%  procedure LDA  
%  procedure LDA  
%  input:    x  N-vector indicating group membership. Only observations
%               with values 1 or 2 are used for estimation, but the
%               linear discriminant function will be applied to all
%               N observations.
%            y  N by p matrix of variables to be used for classification
%       prior1  if 0, prior probabilities are estimated as the relative
%               sample sizes. If between 0 and 1, this will be the
%               prior probability for group 1.
%  output:   b  vector of linear discriminant function coefficients
%            d  p-variate standard distance.
%          ldf  N-vector of linear discriminant function scores, computed
%               for all observations in the data matrix y.
%         xhat  predicted group membership for all N observations
%        post1  posterior probability for group 1, for all N observations


function [b, d, ldf, xhat, post1] = lda(x, y, prior1) ;
y1 = y(find(x==1), :) ;
[N1 p] = size(y1) ; ybar1 = mean(y1)' ; S1 = cov(y1) ;
y2 = y(find(x==2), :) ;
[N2 p] = size(y2) ; ybar2 = mean(y2)' ; S2 = cov(y2) ;
S = ((N1 - 1) * S1  +  (N2 - 1) * S2) / (N1 + N2 - 2) ;
b = inv(S) * (ybar1 - ybar2) ;
d = sqrt( b' * (ybar1 - ybar2) ) ;
ldf = y * b ;
if prior1 == 0 ; prior1 = N1 / (N1 + N2) ; end ;
z = log(prior1/(1-prior1)) - b' * (ybar1 + ybar2)/2  + ldf ;
post1 = exp(z) ./ (1 + exp(z)) ;
post2 = 1 - post1 ; post = [post1 post2] ;
[sortpost, sortind] = sort(post') ;
xhat = sortind(2, :)' ;
%  end of procedure LDA  
%  end of procedure LDA  



