%% Task2.m
% University of Birmingham
% Machine Learning - Computer Based Test 1
% Author: Andreas Farrenkopf
clear all;close all;

%% Load data and extract the mens 100m data
% Change directory/path for use on other computers
load /Users/Andi/Studium/Birmingham/LH_Machine_Learning/MATLAB/data/olympics

% The olypmic years are stored in column 1 and the respective winning times
% in column 2
x = male100(:,1);
t = male100(:,2);

% Rescale x (See also First Course in Machine Learning)
% Further calculations and plots will be more intuitive with the rescale
x = x - x(1);
x = x./4;

% Plot the data points
figure(1);hold on
plot(x,t,'bo','markersize',10);
% X axis does not show year anymore due to rescale
xlabel('Number of olympic game');
ylabel('Time');

%% 1st Order Model
% Compute values from one year before first olympic until one year after
% last olympic games
% Compute value every 0.01 units
% (See also First Course in Machine Learning)
plotx = [x(1)-1:0.01:x(end)+1]';
X = [];
plotX = [];
% Build up data so it has respective order
for k = 0:1
    X = [X x.^k];
    plotX = [plotX plotx.^k];
end

% Compute model paramteres for whole dataset
w = inv(X'*X)*X'*t;

% Plot the linear model
plot(plotx,plotX*w,'r','linewidth',2)

%% Cross Validation for 1st order model
% The cross validation is to find the best lambda for the 1st order model

% Define parameters needed for cross validation
N = length(x); %27
K = 5; %K-fold CV
% Define the range for lambda between 0 and 0.5
lamb = [0:0.01:0.5];

% Call function cross_validation
cv_loss_1 = cross_validation(lamb,K,N,X,t);

% Plot the average error over lambda
figure(2);
plot(lamb,mean(cv_loss_1,1),'linewidth',2)
xlabel('Lambda');
ylabel('Prediction Error');
title('1st order model','fontsize',15)

%% 4th Order Model
% For explanations on code segments and references, please see section '1st Order Model'
plotx = [x(1)-1:0.01:x(end)+1]';
X = [];
plotX = [];
for k = 0:4
    X = [X x.^k];
    plotX = [plotX plotx.^k];
end
w = inv(X'*X)*X'*t;
figure(1);hold on
plot(plotx,plotX*w,'y','linewidth',2)
legend('Data Points','1st Order','4th Order')

%% Cross Validation for 4th order model
% The cross validation is to find the best lambda for the 4th order model
% For explanations on code segments and references, please see section '1st Order Model'
N = length(x); %27
K = 5;
lamb = [0:0.01:0.5];
cv_loss_4 = cross_validation(lamb,K,N,X,t);
figure(3);
plot(lamb,mean(cv_loss_4,1),'linewidth',2)
xlabel('Lambda');
ylabel('Prediction Error');
title('4th order model','fontsize',15)

%% Function definition for cross validation
% The function performs a K fold cross validation based on the input
% parameters
function [cv_loss] = cross_validation(lamb, K, N, X, t)
sizes = repmat(floor(N/K),1,K);
sizes(end) = sizes(end) + N - sum(sizes);
csizes = [0 cumsum(sizes)]; % Sum of all fold sizes must equal 27 (number of data points N)

% Iterate through the range of lambda
for l = 1:length(lamb)
    lambda = lamb(l);
    % Use K fold cross validation to determine lambda with the lowest
    % prediction error on test data
    for fold = 1:K
        % Code amended from 'First Course in Machine Learning'
        % Partition the data into K folds of equal size
        % foldX contains the data for just one fold (To validate the model)
        % trainX contains all other data (To train the model)

        foldX = X(csizes(fold)+1:csizes(fold+1),:);
        trainX = X;
        trainX(csizes(fold)+1:csizes(fold+1),:) = [];
        foldt = t(csizes(fold)+1:csizes(fold+1));
        traint = t;
        traint(csizes(fold)+1:csizes(fold+1)) = [];

        % Compute model parameters with formulae of regularized least
        % squares approach
        w = inv(trainX'*trainX + size(trainX,1)*lambda*eye(size(trainX,2)))*trainX'*traint;

        % Compute prediction values and respective error
        fold_pred = foldX*w;
        cv_loss(fold,l) = mean((fold_pred-foldt).^2);
    end
end
end