% logistic_regression_demo.m
% Generation of a small dataset with two classes with Gaussian distribution
% Compare the results of using mnrfit and those obtained in theory
% Dr. David Pan, UAH

close all

N = 100;

% Generate data entries for Class 1
m1 = [3, 3]';   % Mean vector
cov1 = [2 1; 1 2]; % Covariance matrix
rng default  % For reproducibility
r1 = mvnrnd(m1,cov1,N);

data_C1 = zeros(N, 2);
data_C1 = r1;
label_C1 = ones(N, 1);

% Generate data entries for Class 2
m2 = [9, 9]';
cov2 = [2,1; 1,2];
rng default  % For reproducibility
r2 = mvnrnd(m2,cov2,N);

data_C2 = zeros(N, 2);
data_C2 = r2;
label_C2 = 2*ones(N, 1);

% Combine data of two classes
data = vertcat (data_C1, data_C2);
label = vertcat (label_C1, label_C2);

% Export the dataset to a CSV file (to be imported by sklearn)
data_label = horzcat(data, label);
writematrix(data_label, 'logistic_regress.csv');

gscatter(data(:,1),data(:,2),label,'rb','ox')
grid

% logistic regression
B = mnrfit(data, label);

% Decision boundary based on coefficients in B
% B(2)*x1 + B(3)*x2 + B(1) = 0, where B(1) is the intercept
% log(P(C1|X)/P(C2|X)) >= 0, if P(C1|X)>=P(C2|X)
x1 = min(data(:,1)): 0.01: max(data(:,1));
x2 =  -(B(2)*x1 + B(1))/B(3);
hold on;
plot(x1,x2)

% MLE estimates
mu1 = mean(data_C1);
mu2 = mean(data_C2);
Sigma1 = cov(data_C1);
Sigma2 = cov(data_C2);
Sigma = 1/2*(Sigma1 + Sigma2);

% Slopes and intercepts based on theoretical results
inv(Sigma)*(mu1-mu2)'
-0.5*mu1*inv(Sigma)*mu1'+0.5*mu2*inv(Sigma)*mu2'

% Determine the class probability
x = mean(data)
prob = mnrval(B, x)
log(prob(1)/prob(2))
B(2)*x(1)+B(3)*x(2)+B(1)

x = data(9,:)
x = data(111,:)
label(111,:)
x = data(135,:)
label(135,:)
prob = mnrval(B, x)

% Calculate the classifcation error performance of the entire dataset
P = mnrval(B, data);

label_pred = (P(:,1)>=0.5) + 2*(P(:,2)>=0.5); 
find(label ~= label_pred)
length(find(label ~= label_pred))/length(label)
