load upf_cv11_linear.txt;   % nacti soubor (radky  a sloupce -> matice)

x = upf_cv11_linear(:,1);       % prvni sloupec -> sloupcovy vektor x
y = upf_cv11_linear(:,2);       % ...
sigma = upf_cv11_linear(:,3);   % ...
N = length(x);                % promenna N = pocet radku (pocet mereni)

errorbar(x,y,sigma,"vertical","o"); % vykresli chybove usecky a kulate body
hold on;    % pri dalsim kresleni nesmazavej co uz je nakresleno

% ------ Modelova funkce --------------------------------------------------
% Linearni model s M parametry muzeme napsat jako
%   lambda(x) = a_1 f_1(x) + a_2 f_2(x) + ... + a_M f_M(x)
% Predpovedi pro body y_1,y_2...y_N lze napsat jako
%   lambda = A * x
% kde lambda a x jsou vektory dimenze N a A je matice M x N.
% Prvek A_ij je potom f_j(x_i)

% linearni (polynom 1. stupne):
model = "a + b x";  % text do legendy grafu
A = [ones(N,1), x]; % matice A, i-ty radek bude [1, x_i]
M = 2;      % pocet parametru

% polynom 3. stupne
%model = "a + b x + c x^2 + d x^3";  % text do legendy grafu
%A = [ones(N,1), x, x.^2, x.^3]; % i-ty radek je [1, x_i, x_i^2, x_i^3]
%M = 4;      % pocet parametru

% harmonicka funkce:
%model = "a + b sin(x) + c cos(x) + d sin(2x) + e cos(2x)";
%A = [ones(N,1), sin(x), cos(x), sin(2*x), cos(2*x)];
%M = 5;

%---- F I T ----------------------------------------------------------------
V = diag(sigma.^2);    % matice kovariance (. = udelej operaci po slozkach)
invV = diag(1./sigma.^2);   % inverzni matice k diagonalni matici V
AVA = A.' * invV * A;   % A.' = transpose(A) tedy transponovana matice
B = AVA \ A.' * invV;   % vyres system linearnich rovnic A = AVA * B
theta = B * y;          % nafitovane parametry theta

% vypis nafitovane parametry, vykresli modelovou funkci:
disp(model);
disp("fitted parameters:");
disp(theta.');          % theta je sloupcovy vektor, vypis ho do radku
yModel = A * theta;     % predpoved modelu v namerenych bodech x
plot(x,yModel);         % vykresli predpoved jako caru do grafu

covTheta = B * V * B.'; % matice kovariance parametru theta
disp("covariance:");
disp(covTheta);
% spocitej standardni odchylky sigma z kovariance:
sigmaTheta = zeros(1,M);
for i = 1:M
    sigmaTheta(i) = sqrt(covTheta(i,i));
end
disp("sigmas: ");
disp(sigmaTheta);

% spocite korelacni matici, cor_ij = cov_ij / (sigma_i * sigma_j)
corTheta = covTheta;
for i = 1:M
        for j = 1:M
            corTheta(i,j) = covTheta(i,j) / (sigmaTheta(i)*sigmaTheta(j));
        end
end
disp("correlation:");
disp(corTheta);

%--------------------------------------------------------------------------
xlabel("x");
ylabel("y");
title("Příklad lineární regrese")
legend("y_i \pm \sigma_i",model);
hold off;   % obrazek je hotovy, uz ho muzes necim prepsat
