clear all
clc

%Define points for x axis
x_low=0;
x_high=5;
x=linspace(x_low,x_high,100);
N=length(x);

%Define points for y axis
k=3;
q=0;
y=k*x+20*rand(1,N)+q;

%Calculation of coefficients using linear LSM
N=length(x);
k_LSM=(1/(N^2*var(x)))*(N*sum(x.*y)-sum(x)*sum(y))
q_LSM=(1/(N^2*var(x)))*(-sum(x)*sum(x.*y)+sum(x.*x)*sum(y))

% Measure of quality of the correlation process
y_hat=@(x) k_LSM.*x+q_LSM;

%Deviations of the observations from their mean
SSyy=sum((y-mean(y)).^2); %Note: the expression refers to variance, where 1/N term has been omitted since it will be eliminated 

%Deviations of their observations from their predicted values
SSE=sum((y_hat(x)-y).^2); %Note: the expression refers to variance, where 1/N term has been omitted since it will be eliminated

Determination_index_R_squared=(SSyy-SSE)/SSyy

%Respectively
%Determination_index_R_squared=1-SSE/SSyy

%Plotting the results
figure
plot(x,y,'.','DisplayName','Points')
hold on
fplot(@(x) k_LSM.*x+q_LSM,'DisplayName','Least Squares line fit')
xlim([x_low x_high])
%title(strcat('deviation_{mean}=',num2str(deviation_mean),'; deviation_{std}=',num2str(deviation_std)))
legend show