clear all
clc

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

%Define points for y axis (y(x)=A*x^2+B*x+C)
A=3;
B=2;
C=1;
y=(A*x.^2+B*x+C);
y=y+0.8.*y.*rand(1,N) ;

M=[sum(x.^4) sum(x.^3) sum(x.^2);
   sum(x.^3) sum(x.^2) sum(x);
   sum(x.^2) sum(x) N];
k=[sum(x.^2.*y);sum(x.*y);sum(y)]

b=inv(M)*k;

deviation_mean=1/N*sum(y-(b(1)*x.^2+b(2)*x+b(3)));
deviation_std=sqrt(1/N*sum((y-(b(1)*x.^2+b(2)*x+b(3))).^2));

% Measure of quality of the correlation process
y_hat=@(x) b(1)*x.^2+b(2)*x+b(3);

%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')
xlabel('x')
ylabel('y(x)')
hold on
fplot(@(x) b(1)*x.^2+b(2)*x+b(3),'DisplayName','y(x)=ln(y)=Ax^2+Bx+C')
legend show
% title(strcat('deviation_{mean}=',num2str(deviation_mean),'; deviation_{std}=',num2str(deviation_std)))
xlim([x_low x_high])