
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% EE 675 - Fall 2015
%%%% Simple Code to show How to classify simulated data sets 
%%%%
%%%% Things to try:
%%%% 1. Effect of varying group means m1 and m2
%%%% 2. Effect of varying group standard deviations sd1 and sd2
%%%% 3. Effect of varying number of features
%%%% 4. Effect of varying the classifier type / parameters
%%%% 5. Effect of varying number of training and test samples
%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nf= 1;                      % number of features used
Ntrain= 1000;               % number of training samples
Ntest=  1000;               % number of training samples
m1=     100;                % mean of Group 1 samples
m2=     120;                % mean of Group 2 samples
sd1=    10;                 % std dev of Group 1 samples
sd2=    10;                 % std dev of Group 2 samples

% x = Features
% y = Group labels 

%%%%% Generate samples for each group (training and testing samples)
x1_train=   randn(Ntrain,Nf)*sd1 + m1;       % Group 1 training samples
x2_train=   randn(Ntrain,Nf)*sd2 + m2;       % Group 2 training samples
x1_test=    randn(Ntest,Nf)*sd1 + m1;        % Group 1 testing samples
x2_test=    randn(Ntest,Nf)*sd2 + m2;        % Group 2 testing samples

%%%%% Labels for each group
y1_train=   ones(Ntrain,1);                 % Group 1 train label = 1
y1_test=    ones(Ntest,1);                  % Group 1 test label = 1
y2_train=   2*ones(Ntrain,1);               % Group 2 train label = 2
y2_test=    2*ones(Ntest,1);                % Group 2 test label = 2


%%%% combine all training samples/labels each in one array 
x_train= [x1_train ; x2_train];
y_train= [y1_train ; y2_train];

%%%% combine all training samples/labels each in one array 
x_test= [x1_test ; x2_test];
y_test= [y1_test ; y2_test];


%%%% Linear discriminant analysis classifier
[c1,err,P1] = classify(x_test,x_train,y_train,'linear');     % LDA - Linear
c2 = classify(x_test,x_train,y_train,'quadratic');  % QDA - Quadratic
c3 = knnclassify(x_test,x_train,y_train',1);        % KNN K=1
c4 = knnclassify(x_test,x_train,y_train',1);        % KNN K=3
c5 = knnclassify(x_test,x_train,y_train',1);        % KNN K=3
svmStruct = svmtrain(x_train,y_train);              % SVM - linear (default)
c6 = svmclassify(svmStruct,x_test);                 % SVM
svmStruct = svmtrain(x_train,y_train,'kernel_function','quadratic'); % SVM - quadratic
c7 = svmclassify(svmStruct,x_test);                 % SVM
svmStruct = svmtrain(x_train,y_train,'kernel_function','rbf'); % SVM - RBF
c8 = svmclassify(svmStruct,x_test);                 % SVM
svmStruct = svmtrain(x_train,y_train,'kernel_function','mlp'); % SVM - mlp
c9 = svmclassify(svmStruct,x_test);                 % SVM
c10 = classify(x_test,x_train,y_train,'diaglinear');     % Naive Bayes


err1 = sum(y_test~= c1)/Ntest;         % mis-classification rate
err2 = sum(y_test~= c2)/Ntest;         % mis-classification rate
err3 = sum(y_test~= c3)/Ntest;         % mis-classification rate
err4 = sum(y_test~= c4)/Ntest;         % mis-classification rate
err5 = sum(y_test~= c5)/Ntest;         % mis-classification rate
err6 = sum(y_test~= c6)/Ntest;         % mis-classification rate
err7 = sum(y_test~= c7)/Ntest;         % mis-classification rate
err8 = sum(y_test~= c8)/Ntest;         % mis-classification rate
err9 = sum(y_test~= c9)/Ntest;         % mis-classification rate
err10 = sum(y_test~= c10)/Ntest;         % mis-classification rate

Error_Rate= [err1 , err2 , err3 , err4 , err5 , err6 , err7 , err8 , err9, err10]


%%%%%% Sample confusion matrix computation 
%%% check that the total of the off-diagonal elements = error rate
conMat1 = confusionmat(y_test,c1)/Ntest      % confusion matrix for c1
TP= conMat1(1,1);
TN= conMat1(2,2);
FP= conMat1(1,2);
FN= conMat1(2,1);

Sens= TP/(TP+FN)
Spec= TN/(TN+FP)
PPV= TP/(TP+FP)
NPV= TN/(TN+FN)
Acc= (TP+TN)/(TP+TN+FP+FN)
FNrate= (1- Sens)
FPrate= (1- Spec)
PLR= Sens/(1-Spec)
NLR= (1-Sens)/Spec
DPrev= (TP+FN)/(TP+FN+TN+FP) 






figure(1);
[X,Y,T,AUC] = perfcurve(y_test,P1(:,2),2);
plot(X,Y)
xlabel('False positive rate'); ylabel('True positive rate')
title('ROC for Classification')
AUC


