% 
% Written by:
% -- 
% John L. Weatherwax                2007-07-01
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; clc; clear; 

rand('seed',0);
randn('seed',0);

mu1 = [ 1, 1 ].';
mu2 = [ 1.5, 1.5 ].';
sigmaCov = [ 1.01, 0.2; 0.2, 1.01 ];

nFeats = 10000; 

X1 = mvnrnd( mu1, sigmaCov, nFeats ); 
X2 = mvnrnd( mu2, sigmaCov, nFeats ); 

h1 = plot( X1(:,1), X1(:,2), '.b' ); hold on; 
h2 = plot( X2(:,1), X2(:,2), '.r' ); hold on; 
legend( [h1,h2], {'class 1', 'class 2'} ); 

X = [ X1; X2 ]; 
labels = [ ones(nFeats,1); 2*ones(nFeats,1) ]; 

% Part (a): classify each of the vectors in X1 and X2 (directly using the likelihood ratio)
%
lratio = mvnpdf( X, mu1.', sigmaCov ) ./ mvnpdf( X, mu2.', sigmaCov ); 

% If lratio > 1 we are selecting this sample from class #1 otherwise from class #2: 
%
class_decision = lratio > 1;
choosen_class = zeros(2*nFeats,1);
choosen_class(find(class_decision==1)) = 1;
choosen_class(find(class_decision~=1)) = 2;

P_correct = sum(choosen_class == labels)/(2*nFeats);
P_error   = 1 - P_correct;

% Calculate the optimal Bayes error rate (using the results from Problem~2.9): 
%
addpath('../../../Duda_Hart_Stork/Code/Chapter2/ComputerExercises'); 
dm = mahalanobis(mu1,mu2,sigmaCov);
P_B = 1 - normcdf( 0.5 * dm );

fprintf('empirical P_e= %10.6f; analytic Bayes P_e= %10.6f\n',P_error,P_B); 

% Part (b): classify each of the vectors in X1 and X2 (directly using the likelihood ratio)
%

% If lratio > 1/2 we are selecting this sample from class #1 otherwise from class #2: 
%
class_decision = lratio > 1/2;
choosen_class = zeros(2*nFeats,1);
choosen_class(find(class_decision==1)) = 1;
choosen_class(find(class_decision~=1)) = 2;

P_correct = sum(choosen_class == labels)/(2*nFeats);
P_error   = 1 - P_correct;

% extract the expected loss using this classifier: 
%
L_12 = sum( choosen_class(1:nFeats)~=1 );
L_21 = sum( choosen_class(nFeats+1:end)~=2 );
r_hat = ( L_12 + 0.5 * L_21 )/(2*nFeats);

fprintf('P_correct= %10.6f; r_hat= %10.6f\n',P_correct,r_hat);