% function PR_data_analysis

%% exemplary data preparation
load fisheriris
label=zeros(size(species,1),1);
label(ismember(species,'setosa'))=1;
label(ismember(species,'virginica'))=2;
label(ismember(species,'versicolor'))=3;
data=meas;

%% Linear Discriminant Analysis
% Create a default (linear) discriminant analysis classifier:
linclass = ClassificationDiscriminant.fit(data,label);
% 
% Classify the data with average measurements:
meanmeas = mean(data);
meanclass = predict(linclass,meanmeas);


% Quadratic Discriminant Analysis
quadclass = ClassificationDiscriminant.fit(data,label,...
    'discrimType','quadratic');

% Classify an iris with average measurements using the quadratic classifier:
meanclass2 = predict(quadclass,meanmeas)

%% Linear discriminant analysis
% ldaClass = classify(meas(:,1:2),meas(:,1:2),group);
ldaClass = classify(data,data,label);

% numerical eveluation of the obtained results
bad = eq(ldaClass,label);
ldaResubErr = sum(bad)/length(label);
disp(['The classification rate is ' num2str(ldaResubErr)]);

% comparizon of the obtained results - confusion matrix
[ldaResubCM,grpOrder] = confusionmat(label,ldaClass);  %the diagonal elements represent correctly classified observations
%the (i,j) element in the confusion matrix is the number of samples whose known class label is class i and whose predicted class is j

% if we would like to visualize the regions (it is possible for 2 features only)
% The function has separated the plane into regions divided by lines, and assigned different regions to different species. 
% One way to visualize these regions is to create a grid of (x,y) values and apply the classification function to that grid.

% [x,y] = meshgrid(4:.1:8,2:.1:4.5);
[x,y] = meshgrid(1:.1:10,1:.1:10);
x = x(:);
y = y(:);
j = classify([x y],data(:,2:3),label); % you can change the features to visualize the obtained results
gscatter(x,y,j,'grb','sod');

%% Quadratic discriminant analysis.

qdaClass = classify(data,data,label,'quadratic');
bad = ~eq(ldaClass,label);
ldaResubErr = sum(bad)/length(label);
disp(['The classification error is ' num2str(ldaResubErr)]);

% generate 10 disjoint stratified subsets
cp = cvpartition(label,'k',10);
% Estimate the true test error for LDA using 10-fold stratified cross-validation
% In this case you don't have another labeled data set, 
% but you can simulate one by doing cross-validation. 
% A stratified 10-fold cross-validation is a popular choice 
% for estimating the test error on classification algorithms. 
% It randomly divides the training set into 10 disjoint subsets. 
% Each subset has roughly equal size and roughly the same class 
% proportions as in the training set. Remove one subset, 
% train the classification model using the other nine subsets, 
% and use the trained model to classify the removed subset. 
% You could repeat this by removing each of the ten subsets one at a time.

ldaClassFun= @(xtrain,ytrain,xtest)(classify(xtest,xtrain,ytrain));
ldaCVErr  = crossval('mcr',data,label,'predfun', ...
             ldaClassFun,'partition',cp)
         
% Estimate the true test error for QDA using 10-fold stratified cross-validation.
qdaClassFun = @(xtrain,ytrain,xtest)(classify(xtest,xtrain,ytrain,...
              'quadratic'));
qdaCVErr = crossval('mcr',data,label,'predfun',...
           qdaClassFun,'partition',cp)      
  
%% Gaussian Mixture Model

obj = gmdistribution.fit(data,3)
idx = cluster(obj,data(3,:))

%% test the methodology for the whole avaliable database:

ttest=zeros(size(label));
for j=1:length(ttest)
    dtest=data;
    dtest(j,:)=[];
    dlabel=label;
    dlabel(j,:)=[];
    % teach the classifier
    linclass = ClassificationDiscriminant.fit(dtest,dlabel);
    ttest(j) = predict(linclass,data(j,:));
end

% compare the obtained results
bad = ~eq(ttest,label);
err = sum(bad) / length(label)
%% NaiveBayes Classifier

nbGau= NaiveBayes.fit(data, label);
nbGauClass= nbGau.predict(data);
bad = ~eq(nbGauClass,label);
nbGauResubErr = sum(bad) / length(label)
nbGauClassFun = @(xtrain,ytrain,xtest)...
               (predict(NaiveBayes.fit(xtrain,ytrain), xtest));
nbGauCVErr  = crossval('mcr',data,label,...
              'predfun', nbGauClassFun,'partition',cp)
          
% model each variable in each class using a kernel density estimation, which is a more flexible nonparametric technique.
nbKD= NaiveBayes.fit(data, label,'dist','kernel');
nbKDClass= nbKD.predict(data);
bad = ~strcmp(nbKDClass,label);
nbKDResubErr = sum(bad) / length(label)
nbKDClassFun = @(xtrain,ytrain,xtest)...
            (predict(NaiveBayes.fit(xtrain,ytrain,'dist','kernel'),xtest));
nbKDCVErr = crossval('mcr',data,label,...
            'predfun', nbKDClassFun,'partition',cp)
        
%% Classification Tree
ctree = ClassificationTree.fit(data,label);
view(ctree)   
view(ctree,'mode','graph')
% Resubstitution Error of a Classification Tree
resuberror = resubLoss(ctree)
% predict the classification
Ynew = predict(ctree,data);

OptimalTree = ClassificationTree.fit(data,label,'minleaf',40);
view(OptimalTree,'mode','graph')
YnewOpt = predict(ctree,data)

% prunning
% Find the optimal pruning level by minimizing cross-validated loss:
% [~,~,~,bestlevel] = cvLoss(ctree,...
%     'subtrees','all','treesize','min')

[~,~,~,bestlevel] = cvLoss(ctree,'subtrees','all')
tree = prune(ctree,'Level',bestlevel);
view(tree,'mode','graph')

%% kNN
mdl = ClassificationKNN.fit(data,label,'NumNeighbors',3,...
    'NSMethod','exhaustive','Distance','minkowski');
flwrClass = predict(mdl,data(80,:));

% visualization
x = data(:,3:4);
gscatter(x(:,1),x(:,2),label);
set(legend,'location','best');

newpoint = [5 1.45];
line(newpoint(1),newpoint(2),'marker','x','color','k',...
   'markersize',10,'linewidth',2);

% Find the 10 sample points closest to the new point:
[n,d] = knnsearch(x,newpoint,'k',10)
line(x(n,1),x(n,2),'color',[.5 .5 .5],'marker','o',...
    'linestyle','none','markersize',10)

%% Fuzzy C-Means Clustering

[center,U,obj_fcn] = fcm(data,3); 
[ind,kl]=max(U);
kl=kl';

%% PCA

[coeff,score,latent,tsquared,explained]= pca(data);
explained
[center,Uc,obj_fcn] = fcm(score(:,1:2),3); 
[ind,klc]=max(Uc);
klc=klc';

