Hi Doreen,

Naja mal gucken, ob mir da noch was sinnvolles einfällt.

Ich habe mal schnell das Beispiel aus der Dissertation von Wicki (Kapitel 18.2.4) in Matlab implementiert. Ich poste hier mal den Quelltext in der Hoffnung, dass es Dir hilft. Ich habe es nicht kommentiert und ggf. sind einige Schritte überflüssig. Ich hoffe dennoch, dass Du es lesen kannst.

Viele Grüße
Micha

clear variables;
close all;
format long g;
clc;
 
set(groot,'defaultAxesTickLabelInterpreter','latex');  
set(groot,'defaulttextinterpreter','latex');
set(groot,'defaultLegendInterpreter','latex');
 
% Ausgleichung ohne Fehler in den Beobachtungen
A = [
    -1  1  0  0
     0 -1  1  0
     0 -1  0  1
     0  0 -1  1
     1  0  0  0
     0  0  0  1
     0  0  1  0
     0  1  0  0
    -1  0  1  0
];
 
P = diag(1./0.001^2.*[0.1276 0.1372 0.0772 0.1041 0.0693 0.1041 0.0918 0.1372 0.0865]);
 
y = [
     32.059
     -6.556
     26.170
     32.726
    -27.809
     30.419
     -2.317
      4.246
     25.496
];
 
N = A'*P*A;
 
Qx = inv(N);
x = Qx * A'*P*y;
 
v =  A*x - y;
sigma_v = sqrt(diag(inv(P) - A*Qx*A'));
 
w = v./sigma_v;
 
% Beobachtungsfehler in 1 und 7
c = 3.5;
y([1 7]) = y([1 7]) + [+0.1; -0.1];
 
x = Qx * A'*P*y;
v = A*x - y;
 
Q_v = inv(P) - A*Qx*A';
R   = Q_v * P; 
sigma_v = sqrt(diag(Q_v));
r = diag(R);
 
w = v./sigma_v;
k = c*sigma_v;
 
while abs(max(abs(w)) - c) > 0.01
    % groesste Teststatistik
    [~, idx] = max(abs(w));
 
    d = zeros(size(v));
    d(idx) = v(idx) - sign(v(idx)).*k(idx);
    y = y - d./-r;
 
    x = Qx * A'*P*y;
    v = A*x - y;
 
    w = v./sigma_v;
end
 
disp(x);
disp(w);

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences