H_R  = 1.00 ;
h_ob = [ 1.04 2.05 3.03 ]' ;
A    = [ 1 0 ; -1 1 ; 0 1 ]; 
m    = [ .02 ;  .02 ; .04 ];
P = inv(diag(m.^2));
L    = [ -H_R ; 0 ; -H_R ] - h_ob;
X = -inv(A'*P*A) * (A'*P*L);
V = A*X + L;

h_w = h_ob + V;
sigma0 = sqrt(V'*P*V/ (size(A)(1) - size(A)(2)))

% część nowa, ocena dokładności
Cx = sigma0**2 * (A'*P*A)^-1
mH = sqrt(diag(Cx))

rho = Cx(2,1) / prod(mH)

Cv = sigma0^2 * (inv(P) - A*inv(A'*P*A)*A')
mv = sqrt(diag(Cv))

Ch = sigma0**2 * (A*inv(A'*P*A)*A')
mh = sqrt(diag(Ch))

D = [-1 1]
Cu = D * Cx * D'
mu = sqrt(Cu)

D = [ 0 2 ; 2  -1]
Cu = D * Cx * D'

D = [1 1 1 ; 1 -1 0 ]
Cu= D*Ch*D'
mu = sqrt(diag(Cu))
