%% Lab 9
% Based on the Classification homework problem from 2017

clear all
close all
load('fisherData');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Part 1 - Prototypes (equiv. to assuming circular covariance/zero correlation)

ndata1 = size(data1,1);     % number of examplars for class 1
ndata2 = size(data2,1);     % class 2
mu1 =                       % compute class means
mu2 = 
vproto =                    % discriminant vector between the two means
mumean =                    % assume criterion halfway between the means
thresh =                    % compute threshold for discriminant
vorthog =                   % compute vector of length 1 orthogonal to discriminant

figure(1);
subplot(2,2,1);
plot(data1(:,1),data1(:,2),'r.');
hold on;
plot(mu1(1), mu1(2),'rs','MarkerSize',12);
plot(data2(:,1),data2(:,2),'g.');
plot(mu2(1), mu2(2),'gs','MarkerSize',12);
plot(...)                   % plot discriminant boundary
title('prototype discriminant');
axis equal

% run the discriminant on the data on which it was trained
% compute "hits","false alarms", etc., percent correct, d'
% Why is d' a silly statistic?

%% Part 2 - Fisher linear discriminant
avgcov =                    % compute pooled/averaged covariance
vfisher =                   % compute Fisher discriminant vector
thresh = 
vorthog = 

subplot(2,2,2);
plot(data1(:,1),data1(:,2),'r.');
hold on;
plot(mu1(1), mu1(2),'rs','MarkerSize',12);
plot(data2(:,1),data2(:,2),'g.');
plot(mu2(1), mu2(2),'gs','MarkerSize',12);
plot(...)                   % plot discriminant boundary
title('Fisher discriminant');
axis equal

% run the discriminant on the data on which it was trained as above

%% Part 3 - QDA log likelihoods

[xgrid,ygrid] = meshgrid(...);      % compute an xy grid of sample locations
ll1 =                               % compute log likelihoods on grid for category 1
ll2 =                               % same for category 2
subplot(2,2,3);
surf(xgrid,ygrid,reshape(ll1,size(xgrid)));
hold on;
surf(xgrid,ygrid,reshape(ll2,size(xgrid)));
title('log likelihoods');

% run the QDA on the data

ll1data1 =                          % log likelihood for category 1 data vs. category 1 Gaussian
ll2data1 =                          % same for category 2's Gaussian
llrdata1 = ll2data1 - ll1data1;     % log likelihood ratio
inds1data1 = (llrdata1 < 0);        % indices that will be misclassified as category 2
inds2data1 = (llrdata1 > 0);        % or correctly as category 1, I suppose equal to 0 needs to be assigned
% do the same for category 2 data

subplot(2,2,4);
plot(data1(inds1data1,1),data1(inds1data1,2),'r.');
hold on;
plot(data1(inds2data1,1),data1(inds2data1,2),'rx');
plot(mu1(1), mu1(2),'rs','MarkerSize',12);
plot(data2(inds1data2,1),data2(inds1data2,2),'gx');
plot(data2(inds2data2,1),data2(inds2data2,2),'g.');
plot(mu2(1), mu2(2),'gs','MarkerSize',12);
axis equal

xrg = get(gca,'xlim');  yrg = get(gca, 'ylim');
npts = 200;
[X,Y] = meshgrid(linspace(xrg(1),xrg(2),npts),linspace(yrg(1),yrg(2),npts));
ll1 =                               % same for this new grid
ll2 = 
contour(X,Y,reshape(ll1-ll2,npts,npts), [0, 0],'b');
set(gca,'xlim',xrg); set(gca,'ylim',yrg);
title('QDA');
axis equal;
hold off;

% compute same stats

% Other things you can try:
%
% Cross validation (leave-one-out)
%
% Construct other example datasets sampled from pairs of means and covariances. 
% Compare the discrimination boundaries for each discriminant above to the theoretical 
% QDA boundary from the theoretical underlying distribution
