% program to generate data from two classes and to implement a classifier
% based on assuming Gaussian class-conditional distribution
clear all;
% data generation (use 0.5, 1, and 2 in datagen)
nTr = 200; 
nTe = 400;
[dX dY] = cl_data(nTr,[],0);
ii_1  = find(dY==1);
data1 = dX(ii_1,:);
data2 = dX(setdiff([1:nTr],ii_1),:);
N_1   = size(data1,1);
N_2   = size(data2,1);
pi1 =N_1/nTr; pi2 = N_2/nTr;

% means (mju) and covariances (sig) of the datasets:
mu1 = sum(data1,1)/N_1;  mu2 = sum(data2,1)/N_2;

d = data1 - repmat(mu1,[N_1 1]);
sig1 = d'*d/N_1;  sig1inv = sig1^(-1);
d = data2 - repmat(mu2,[N_2 1]);
sig2 = d'*d/N_2;  sig2inv = sig2^(-1);

% variables used for generating graphs
minx = -0.0; maxx = 5.0; miny = -0; maxy = 6.5;
delta = 0.05;
[gridx,gridy] = meshgrid(minx:delta:maxx,miny:delta:maxy);

% computing the class-conditional probabilities
points = [gridx(:) gridy(:)];
p1     = points - repmat(mu1,[length(points) 1]);
temp   = (p1*sig1inv).*p1;
temp   = exp( -0.5 * sum(temp,2));
c1     = pi1*reshape(temp,size(gridx))/sqrt(det(2*pi*sig1));

p2     = points - repmat(mu2,[length(points) 1]);
temp   = (p2*sig2inv).*p2;
temp   = exp( -0.5 * sum(temp,2));
c2     = pi2*reshape(temp,size(gridx))/sqrt(det(2*pi*sig2));

% COMPUTING posterior probabilities
postc1 = c1./ (c1+c2);

% graphic visualisation
pointSize = 20; figure(1); clf;
colormap gray; subplot 221; hold off;
imagesc(gridx(1,:),gridy(:,1), c1); axis xy;
hold on;
scatter(data1(:,1),data1(:,2),pointSize,[0.7 0.7 0],'d');
title('Class-cond. prob. C_1');
subplot 222; hold off;
imagesc(gridx(1,:),gridy(:,1), c2); axis xy;
hold on;
scatter(data2(:,1),data2(:,2),pointSize,[1 0 0],'o');
title('Class-cond. prob. C_2');

% generating graphics for posterior probability
subplot 223; hold off;
imagesc(gridx(1,:),gridy(:,1),postc1); axis xy;
title('Posterior probability'); 
subplot 224; hold off
[c,h] = contour(gridx(1,:),gridy(:,1),postc1,[0.5 0.5]);
set(h,'Linewidth',3,'Edgecolor',[0.2 0.4 0.1]); hold on
scatter(data1(:,1),data1(:,2),pointSize,[0 0 0],'d');
scatter(data2(:,1),data2(:,2),pointSize,[.5 0 0],'o');
title('Classification boundary');

% generating nTe test data
[dX, dY] = cl_data(nTe,[],0);
ii_1  = find(dY==1);
data1 = dX(ii_1,:);
data2 = dX(setdiff([1:nTe],ii_1),:);
N_1   = size(data1,1);
N_2   = size(data2,1);

% performing classification
log_sig1 = log(det(sig1))-log(pi1);
log_sig2 = log(det(sig2))-log(pi2);

% correct/incorrect classification for data1
p_1 = data1 - repmat(mu1,[length(data1) 1]);
p_1 = sum(((p_1*sig1inv).*p_1),2);

p_2 = data1 - repmat(mu2,[length(data1) 1]);
p_2 = sum(((p_2*sig2inv).*p_2),2);

class(1,1) = length( find(-p_1+p_2 > log_sig1-log_sig2) );
class(1,2) = N_1 - class(1,1);

% SIMILARLY for the second dataset
p_1 = data2 - repmat(mu2,[length(data2) 1]);
p_1 = sum(((p_1*sig2inv).*p_1),2);

p_2 = data2 - repmat(mu1,[length(data2) 1]);
p_2 = sum(((p_2*sig1inv).*p_2),2);

class(2,2) = length( find( -p_1+p_2 > log_sig2-log_sig1));
class(2,1) = N_2 - class(2,2);

% printing out the correct classifications and the error rates for the
% two classes:
class
