clear all;
clc;
%%  Vyrovnávání dat
A = [1   -1   -1   0    0    0;
     0    1    0  -1    0    0;
     0    0    1   0   -1    0;
	 0    0    0   1    1   -1];
 
dat = [101.91 68.45 34.65 64.20 36.44 98.88];
x = dat';
 
V = diag([1 1 1 1 1 1]);
 
W = eye(length(x)) - V*A'*inv(A*V*A')*A;
x_rec = W*x;
V_rec = W*V*W';
std_rec = sqrt(diag(V_rec));

%%  Detekce hrubých chyb

%vektor reziduí
r = A*x;
V_r = A*V*A';
std_V_r = sqrt(diag(V_r));

%   Global Test (GT)
gama = r'*inv(V_r)*r;               %Hodnota statistiky GT je rovna hodnotě minima objektivní funkce u vyrovnávání dat.

%   Constraint nebo Nodal Test (NT)
z_CT = abs(r)./sqrt(diag(V_r));          % z_r(i) = |r(i)|/sqrt(V_r(ii))
%gausswin

%   Measurement Test (MT)
v_MT = V*A'*inv(V_r)*r;             % v = y - y_hat
%V_v = V*A'*inv(V_r)*V;              %z kapitoly 7 z Data Reconciliation and Gross Error detection
V_v = V*A'*inv(V_r)*A;              %!!! upravil jsem, ale vychází moc velká čísla - patně
z_MT = abs(v_MT)./sqrt(diag(V_v));
d = inv(V)*v_MT;
%V_d = A'*inv(A*V*A')*A';       %z kapitoly 7 z Data Reconciliation and Gross Error detection
V_d = A'*inv(A*V*A')*A;         % !!! kdybych to ale upravil takto
z_MT_d = abs(d)./sqrt(diag(V_d));

%   Pro alfa = 5% je ''kritická'' hodnota 1,96 - musím to přepsat podle správné terminologie!


