% LFD 2005. Tutorial Week 4 Sample Answer: Gaussian Classifier
% Tutor: Timothy Hospedales

load('week4.mat');
y=xtest;
figure(1); clf;
plot3(x1(:, 1), x1(:, 2), x1(:, 3), 'b.', x2(:, 1), x2(:, 2), x2(:, 3), 'r.'); grid on; hold on; title('Training Data'); legend('Train c1', 'Train c2');

%Prior class probability
n1 = size(x1, 1); 
n2 = size(x2, 1);
nt = size(y, 1);
pr1 = n1/(n1+n2);
pr2 = n2/(n1+n2);

%ML estimates for means
x1m = mean(x1, 1)';
x2m = mean(x2, 1)';

%Use three attributes
%ML estimate for covariance/precision
c1 = cov(x1); p1 = inv(c1);
c2 = cov(x2); p2 = inv(c2);

%Posterior class probability: exp(-0.5(x-m)*(c^-1)*(x-m))/(2*pi*det(c1))
pc1_y = -0.5*sum((y - repmat(x1m', nt,1)) .* (p1 * (y - repmat(x1m', nt,1))')',2) -0.5*log(det(c1)) + log(pr1);
pc2_y = -0.5*sum((y - repmat(x2m', nt,1)) .* (p2 * (y - repmat(x2m', nt,1))')',2) -0.5*log(det(c2)) + log(pr2);

%Make a decision based on posterior prob
cl1 = find(pc1_y > pc2_y);
cl2 = find(pc1_y < pc2_y);
figure(2); clf; plot3(y(cl1, 1), y(cl1, 2), y(cl1, 3), 'bx', y(cl2, 1), y(cl2, 2), y(cl2, 3), 'rx'); 
plot3(x1(:, 1), x1(:, 2), x1(:, 3), 'b.', x2(:, 1), x2(:, 2), x2(:, 3), 'r.', y(cl1, 1), y(cl1, 2), y(cl1, 3), 'bx', y(cl2, 1), y(cl2, 2), y(cl2, 3), 'rx');grid on;
legend('Train c1', 'Train c2', 'Test c1', 'Test c2'); title('Test Data: 3 Attribs');

%Use two attributes
%ML estimate for covariance/precision
c1 = cov(x1(:,1:2)); p1 = inv(c1);
c2 = cov(x2(:,1:2)); p2 = inv(c2);

%Posterior class probability: exp(-0.5(x-m)*(c^-1)*(x-m))/(2*pi*det(c1))
pc1_y = -0.5*sum((y(:,1:2) - repmat(x1m(1:2)', nt,1)) .* (p1 * (y(:,1:2) - repmat(x1m(1:2)', nt,1))')',2) -0.5*log(det(c1)) + log(pr1);
pc2_y = -0.5*sum((y(:,1:2) - repmat(x2m(1:2)', nt,1)) .* (p2 * (y(:,1:2) - repmat(x2m(1:2)', nt,1))')',2) -0.5*log(det(c2)) + log(pr2);

%Make a decision based on posterior prob
cl1 = find(pc1_y > pc2_y);
cl2 = find(pc1_y < pc2_y);
figure(3); clf; plot(y(cl1, 1), y(cl1, 2), 'bx', y(cl2, 1), y(cl2, 2), 'rx'); 
plot(x1(:, 1), x1(:, 2), 'b.', x2(:, 1), x2(:, 2), 'r.', y(cl1, 1), y(cl1, 2),'bx', y(cl2, 1), y(cl2, 2), 'rx');
legend('Train c1', 'Train c2', 'Test c1', 'Test c2'); title('Test Data: 2 Attribs');

% Or, for non-vectorized class probability
% for i = 1 : nt
%     pc1_y(i) = -(y(i,:)'-x1m)' * p1 * (y(i,:)'-x1m) - 0.5*log(det(c1)) + log(pr1);
%     pc2_y(i) = -(y(i,:)'-x2m)' * p2 * (y(i,:)'-x2m) - 0.5*log(det(c2)) + log(pr2);
% end