% Perform logistic regression for a 2-class dataset. Learn regression params with built-in MATLAB/Octave unconstrained
% optimization routine, fminunc.
% Overview:
% Take in a csv file (without header)
% Split into a training and test set
% Minimize cost function with fminunc
% Compute performance metrics:
%	prediction accuracy on training and test sets
%	confusion matrix, specificity, sensitivity 
%	misclassification error

clc;clear all;close all;

%Input must contain feature columns followed by dependent variable column at end
data = load('class_function_01.txt');

%percentage of data to use for training
train_frac = .75;

%threshold for classifying hypothesis output
thresh = 0.5;

%extract columns to use
X = data(:,1:end-1);
y = data(:,end);

%split into training and test sets:
test_rows = round(size(X,1)*(1-train_frac)); %number of rows to use in test set
X_test = X(1:test_rows,:);
y_test = y(1:test_rows,:);%this is the test set
X_train = X(test_rows+1:end,:);
y_train = y(test_rows+1:end,:);%this is the training set

%Add intercept term to X
X_train = [ones(size(X_train,1), 1) X_train];
X_test = [ones(size(X_test,1), 1) X_test];

% Initialize fitting parameters
initial_thetas = zeros(size(X_train,2), 1);

%  Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 100);

%  Run fminunc to obtain the optimal theta
%  fminunc finds a local minimum of a function of several variables.
[thetas, cost] = fminunc(@(t)(costFunction(t, X_train, y_train)), initial_thetas, options);

% Print theta to screen
fprintf('Cost at theta found by fminunc: %f\n', cost);
fprintf('theta: \n');
fprintf(' %f \n', thetas);

fprintf('\n*********************\n');
%prediction accuracy
p_train = double(logsig(X_train*thetas) >= thresh);
fprintf('Training set accuracy: %f\n', mean(double(p_train == y_train)) * 100);

p_test = double(logsig(X_test*thetas) >= thresh);
fprintf('Test set accuracy: %f\n', mean(double(p_test == y_test)) * 100);

%confusion matrix, sensitivity, specificity, misclassification error
cm = confMatrix(y_test,p_test);
sens = cm(1,1) / (cm(1,1) + cm(1,2)); %ability to identify positive class
spec = cm(2,2) / (cm(2,2) + cm(2,1)); %ability to identify negative class

testError = misclassError(y_test,logsig(X_test*thetas),thresh); %0/1 misclassification error on test set

fprintf('\n*********************\n');
fprintf('Confusion Matrix:\n');disp(cm)
fprintf('Sensitivity: %g\n',sens);
fprintf('Specificity: %g\n',spec);
fprintf('Misclassification error on test set: %g\n',testError);
