function  bayes_classifier(feature1,feature2,p1,p2);
 
%-----------------------------------------------------
% bayes_classifier(class1,class2,p1,p2);
% BAYES_CLASSIFIER function calculates the discriminant functions for
% two classes.
% INPUT variables:
% - feature1 - feature vector of the size (m,2) for the first class (m is the number of samples in the 1st class)
% - feature2 - feature vector of the size (n,2) for the second class (n is the number of samples in the 2nd class)
% - p1 - prior probability for the first class
% - p2 - prior probability for the second class
% OUTPUT:
% Discriminant functions for two classes: 
% g1(x) = A1*x1^2 + B1*x1*x2 + C1*x2^2 + D1*x1 + E1*x2 + F1
% g2(x) = A2*x1^2 + B2*x1*x2 + C2*x2^2 + D2*x1 + E2*x2 + F2

% --------------------------------------------------
% Evgeny Krestyannikov
% krestyan@cs.tut.fi
% Institute of Signal Processing
% Room TE 313

syms x1
syms x2

cov1=cov(feature1);
cov2=cov(feature2);
mu1=mean(feature1);
mu2=mean(feature2);

detcov1=det(cov1);
detcov2=det(cov2);
invcov1=inv(cov1);
invcov2=inv(cov2);

W1=-0.5*invcov1;
w1=invcov1*mu1';

W2=-0.5*invcov2;
w2=invcov2*mu2';


omega1=-0.5*(mu1)*invcov1*mu1'-0.5*log(detcov1)+log(p1);
omega2=-0.5*(mu2)*invcov2*mu2'-0.5*log(detcov2)+log(p2);

omega=omega1-omega2;

X=[sym(x1); sym(x2)];
X_t=[sym(x1) sym(x2)];

digits(2);
a=X_t*(sym(W1,'d')*X);
a=simplify(a);

b=(sym(w1','d')*X);

c=X_t*(sym(W2,'d')*X);
c=simplify(c);

d=sym(w2','d')*X;

eq1=a+b+sym(omega1,'d');
eq2=c+d+sym(omega2,'d');


disp(['Discriminant function for the first class:']);
disp(eq1);
disp(['Discriminant function for the second class:']);
disp(eq2);