Merge remote-tracking branch 'origin/main'

Parameterschaetzung.py

@@ -0,0 +1,26 @@
+import sympy as sp
+
+#für Varianzkomponentenschätzung
+MAX_ITER = 10
+TOL = 1e-3  # 0.1%.
+
+for loop in range(MAX_ITER):
+
+    Q_ll, P = modell.aufstellen_Qll_P()
+
+    N = A.T * P * A
+    n_vec = A.T * P * l
+    dx = N.LUsolve(n_vec)
+
+    v = l - A * dx
+
+    sigma_hat = modell.varianzkomponenten(v, A)
+
+    print(f"Iteration {loop+1}, σ̂² Gruppen:", sigma_hat)
+
+    # Prüfen: ist jede Komponente ≈ 1?
+    if all(abs(val - 1) < TOL for val in sigma_hat.values()):
+        print("Konvergenz erreicht ✔")
+        break
+
+    modell.update_sigma(sigma_hat)

Stochastisches_Modell.py

@@ -3,7 +3,7 @@ from dataclasses import dataclass, field
 from typing import Dict, Tuple
 
 @dataclass
-class StochastischesModellApriori:
+class StochastischesModell:
     sigma_beob: Iterable[float]           #σ der einzelnen Beobachtung
     group_beob: Iterable[int]             #Gruppenzugehörigkeit jeder Beobachtung (Distanz, Richtung, GNSS, Nivellement,...,)
     sigma0_groups: Dict[int, float] = field(default_factory=dict)  #σ0² für jede Gruppe
@@ -32,7 +32,7 @@ class StochastischesModellApriori:
         Q_ll = sp.zeros(n, n)
         P = sp.zeros(n, n)
 
-        for i in range(n):
+        for i in range(self.n):
             sigma_i = self.sigma_beob[i, 0]         #σ-Wert der i-ten Beobachtung holen
             g = int(self.group_beob[i, 0])          #Gruppenzugehörigkeit der Beobachtung bestimmen
             sigma0_sq = self.sigma0_groups[g]       #Den Varianzfaktor der Gruppe holen
@@ -43,70 +43,44 @@ class StochastischesModellApriori:
 
 
     @staticmethod
-    def redundanz_pro_beobachtung(A: sp.Matrix, P: sp.Matrix) -> sp.Matrix:
-        n_beob = P.rows                             #Anzahl der Beobachtungen (Zeilen in P)
-        n_param = A.cols                            #Anzahl der Unbekannten (Spalten in A)
-
-        sqrtP = sp.zeros(n_beob, n_beob)      #Wurzel von P (der Diagonale)
-        for i in range(n_beob):
+    def redundanz_pro_beobachtung(A, P):
+        n = P.rows
+        sqrtP = sp.zeros(n, n)
+        for i in range(n):
             sqrtP[i, i] = sp.sqrt(P[i, i])
 
         A_tilde = sqrtP * A
-
         M = (A_tilde.T * A_tilde).inv()
 
-        r_vec = sp.zeros(n_beob, 1)
-
-        for i in range(n_beob):
-            a_i = A_tilde.row(i)          # 1 × n_param
-            a_i_row = sp.Matrix([a_i])    # explizit 1×n-Matrix
-            r_i = 1 - (a_i_row * M * a_i_row.T)[0, 0]
-            r_vec[i, 0] = r_i
-
-        return r_vec
+        r = sp.zeros(n, 1)
+        for i in range(n):
+            a_i = sp.Matrix([A_tilde.row(i)])
+            r[i] = 1 - (a_i * M * a_i.T)[0]
+        return r
 
 
-    def varianzkomponentenschaetzung(
-        self,
-        v: sp.Matrix,   # Residuenvektor (n × 1)
-        A: sp.Matrix,   # Designmatrix
-    ) -> Dict[int, float]:
+    def varianzkomponenten(self, v, A) -> Dict[int, float]:
+        _, P = self.aufstellen_Qll_P()
+        r_obs = self.redundanz_pro_beobachtung(A, P)
+        gruppen = sorted(set(int(g) for g in self.group_beob))
 
-        if v.rows != self.n_beob:
-            raise ValueError("Länge von v passt nicht zur Anzahl Beobachtungen im Modell.")
+        sigma_hat = {}
 
-        # Aktuelle Gewichte
-        Q_ll, P = self.aufstellen_Qll_P()
+        for g in gruppen:
+            idx = [i for i in range(self.n) if int(self.group_beob[i]) == g]
 
-        # Redundanzzahlen pro Beobachtung
-        r_vec = self.redundanz_pro_beobachtung(A, P)
+            v_i = sp.Matrix([v[i] for i in idx])
 
-        new_sigma0_sq: Dict[int, float] = {}
+            P_i = sp.zeros(len(idx))
+            for k, j in enumerate(idx):
+                P_i[k, k] = P[j, j]
 
-        # Für jede Gruppe j:
-        unique_groups = sorted({int(g) for g in self.group_beob})
+            r_g = sum(r_obs[j] for j in idx)
 
-        for g in unique_groups:
-            # Indizes der Beobachtungen in dieser Gruppe
-            idx = [i for i in range(self.n_beob) if int(self.group_beob[i, 0]) == g]
-            if not idx:
-                continue
+            sigma_hat[g] = float((v_i.T * P_i * v_i)[0] / r_g)
+        return sigma_hat
 
-            # v_j, P_j, r_j extrahieren
-            v_j = sp.Matrix([v[i, 0] for i in idx])               # (m_j × 1)
-            P_j = sp.zeros(len(idx), len(idx))
-            r_j = 0
-            for ii, i in enumerate(idx):
-                P_j[ii, ii] = P[i, i]
-                r_j += r_vec[i, 0]
 
-            # σ̂_j^2 = (v_jᵀ P_j v_j) / r_j
-            sigma_hat_j_sq = (v_j.T * P_j * v_j)[0, 0] / r_j
-
-            # als float rausgeben, kann man aber auch symbolisch lassen
-            new_sigma0_sq[g] = float(sigma_hat_j_sq)
-        return new_sigma0_sq
-
-    def update_sigma0(self, new_sigma0_sq: Dict[int, float]) -> None:
-        for g, val in new_sigma0_sq.items():
-            self.sigma0_groups[int(g)] = float(val)
+    def update_sigma(self, sigma_hat_dict):
+        for g, val in sigma_hat_dict.items():
+            self.sigma0_groups[g] = float(val)