%% Statistics and Probability Applications solved with Matlab
clear all
close all

%% load data from spreadsheet
filename = 'Assignment_Data.xlsx';
sheet = 1;
[nums, txt] = xlsread(filename,sheet,'','basic');% make sure the working directory contains the excel file

% seperate 'number of service call requests' from 'number of days'
NoOfRequests = nums(1,:);
NoOfDays = nums(2,:);


%% group number of requests (range=3)
GroupRange = 3;
GroupedRequests = [];% create empty cell array (3 rows: boundaries; grouped no. of days; mid class values)

Boundary1 = min(NoOfRequests);% set first lower boundary
MaxReached = false;
i = 1;% set first index
while ~MaxReached
    Boundary2 = Boundary1 + GroupRange;% set upper boundary
    % put lower and upper boundaries in cell array
    GroupedRequests(i,1) = Boundary1;
    GroupedRequests(i,2) = Boundary2;
    
    Boundary1 = Boundary2;% set new lower boundary
    i = i+1;% increase index
    
    if Boundary2>=max(NoOfRequests)% if the maximum value has been reached...
        MaxReached = true;% stop setting boundaries
    end
end

%% define vector GNOD
% loop through each request/no. of days and sum all the days that fall into
% each category
GNOD = zeros(1,size(GroupedRequests,1));% initialise GNOD to 0
for i = 1:length(NoOfDays)% loop through number of days
    for j = 1:size(GroupedRequests,1)% loop through request classes
        if NoOfRequests(i)>=GroupedRequests(j,1) && NoOfRequests(i)<GroupedRequests(j,2)
            GNOD(j) = GNOD(j) + NoOfDays(i);% keep track of days in each category
            break% move to next number of requests if the group for this number has been found
        end
    end
end

%% define vector MCV
%sum boundaries of each class, and divide by 2
MCV = sum(GroupedRequests,2)./2;% sum along dimension 2 (i.e. sum each group of boundaries)

%% Plot histogram and polygon
figure;% open new figure
% As 'N' upper boundary is equal to 'n+1' lower boundary, bin edges are
% defined using the lower boundaries, with the final edge defined by the
% final upper boundary
histogram('BinEdges',[GroupedRequests(:,1)' GroupedRequests(end,2)] ,'BinCounts',GNOD,...
            'FaceColor',[0 0 1])% plot histogram with specified options
hold on% required to plot on top of histogram
plot(MCV,GNOD,'r-','Linewidth',1.5)% plot grouped number of days against mid class value

xticks([GroupedRequests(1,1):1.5:GroupedRequests(end,2)])% lowest lower boundary to highest upper boundary
yticks([0:1:max(GNOD)+1])
xlabel('Number of Requests Classes')
ylabel('Number of Days')
title('Histogram of Polygon distribution of data set')

% find modal value
[~,mode_i] = max(GNOD); % find index of modal value
mode_value = MCV(mode_i);

% plot modal value at the modal position along the x axis at a height of 2.
% horizontal alignment of center allows the text to appear in the centre of
% the modal value.
text(mode_value,2,['Modal Value = ' num2str(mode_value)],...
                    'color','y','FontWeight','bold','FontSize',13,'HorizontalAlignment','center')

%% Plot the cumulative + inverse cumulative distributions
% cumulative
cumulative_x = [GroupedRequests(:,1)' GroupedRequests(end,2)]; % use lower boundaries and then final upper boundary again
cumulative_days = zeros(1,length(cumulative_x)); % set expected size of variable
total_sum = sum(GNOD);% find total number of days to "normalise" values

cumulative_days(1) = 0;% set first value to be zero (up to first lower boundary, sum is 0)
cumulative_days(2:end) = cumsum(GNOD);% populate from second index onwards
norm_cumulative_days = cumulative_days./total_sum;% normalise cumulative sum
% interpolate data
interp_x = cumulative_x(1):0.01:cumulative_x(end);% define new x array to interpolate data
interp_norm_cumulative_days = interp1(cumulative_x,norm_cumulative_days,interp_x,'pchip');
% Using 'pchip' interpolated method to avoid straight lines between points,
% pchip produces smooth curve

inverse_norm_cumulative_days = (total_sum-cumulative_days)./total_sum;% total value minus the cumulative value found previously, then normalised
interp_inverse_norm_cumulative_days = interp1(cumulative_x,inverse_norm_cumulative_days,interp_x,'pchip');
% interpolate using previous new x array

figure;
% plot cumulative
plot(cumulative_x,norm_cumulative_days,'bo')
hold on
plot(interp_x,interp_norm_cumulative_days,...
        'color','b','LineWidth',1.5,'LineStyle','--')
% plot inverse cumulative
plot(cumulative_x,inverse_norm_cumulative_days,'ro')
plot(interp_x,interp_inverse_norm_cumulative_days,...
    'color','r','LineWidth',1.5,'LineStyle','-.')


xlim([cumulative_x(1) cumulative_x(end)])
xticks([[cumulative_x(1):cumulative_x(end)]])
grid minor
grid on
yticks([0:0.1:1])
ylabel('Frequencies')
xlabel('xi')
title('Ogive Distribution Frequencies')


% take difference between inverse cumulative and cumulative data
dif_inverseCumulative_Cumulative = interp_inverse_norm_cumulative_days-interp_norm_cumulative_days;
% Where the sign of the difference changes shows the index after the lines cross
intersection_i = find(diff(sign(dif_inverseCumulative_Cumulative)));
% the sign changes because it shows where the larger of the 2 values becomes the smaller of the 2 values 

% take index BEFORE and AFTER cross and take the mean
median_x = mean([interp_x(intersection_i-1) interp_x(intersection_i)]);

% print line of "#" to seperate results
fprintf('\n###########################################################\n')
% start and end with a new line, to seperate it from other text
fprintf(['\nThe median of the data set is ' num2str(median_x) ' call requests.\n'])


%% export data into excel
% setup cell array to write to excel Sheet 2
Sheet2 = cell(8,length(MCV)+1);
% 8 rows (categories to populate)
% columns are the number of classes +1 for row names

% populate categories (col 1)
Sheet2{1,1} = 'Number of Requests xi';
Sheet2{2,1} = 'Mid Class ci';
Sheet2{3,1} = 'Number of days ni';
Sheet2{4,1} = 'Frequency';
Sheet2{5,1} = 'Cumulated Frequency';
Sheet2{6,1} = 'Inversely Cumulated Frequency';
Sheet2{7,1} = 'ni.ci';
Sheet2{8,1} = 'ni.ci^2';

% populate data for all rows, column by column
for i = 1:length(GroupedRequests)
    Sheet2{1,i+1} = ['[' num2str(GroupedRequests(i,:)) ']'];% "+1" to avoid overwriting category
    Sheet2{2,i+1} = MCV(i);
    Sheet2{3,i+1} = GNOD(i);
    Sheet2{4,i+1} = GNOD(i)/total_sum;
    Sheet2{5,i+1} = norm_cumulative_days(1+i);
    Sheet2{6,i+1} = inverse_norm_cumulative_days(i);
    Sheet2{7,i+1} = GNOD(i)*MCV(i);
    Sheet2{8,i+1} = GNOD(i)*(MCV(i)^2);
end

% write cell array (SHEET2) to excel sheet 2
writecell(Sheet2,filename,'Sheet',2)


%% find mean, variance, and standard deviation of the *grouped* data
% create matrix of all "raw" data
all_data_grouped = zeros(1,sum(GNOD)); % create empty matrix to populate
last_i = 0;
for i = 1:length(GNOD)% loop through classes
   all_data_grouped(last_i+1:last_i+GNOD(i)) = MCV(i);% repeat the MCV the amount of times specified in GNOD 
   last_i = last_i+GNOD(i);% used to make sure the next iteration does not overwrite values
end
 
% ##########################################################
% Note: This is not the mean/variance/etc of the total raw 
% data (individual days). This is the stats for the grouped 
% data, as the question asks. The data uses the Mid Class 
% value to find the statistics.
% ##########################################################
% mean
mean_x = mean(all_data_grouped);
fprintf(['\nThe mean of the grouped data set is ' num2str(mean_x) ' call requests.\n'])

% variance
var_x = var(all_data_grouped);
fprintf(['\nThe variance of the grouped data set is ' num2str(var_x) ' call requests.\n'])

%standard deviation
std_x = std(all_data_grouped);
fprintf(['\nThe standard deviation of the grouped data set is ' num2str(std_x) ' call requests.\n'])

% end with a line of '#' to seperate output
fprintf('\n###########################################################\n')