diff --git a/Parameterschaetzung.py b/Parameterschaetzung.py index e69de29..7b394a7 100644 --- a/Parameterschaetzung.py +++ b/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) \ No newline at end of file diff --git a/Stochastisches_Modell.py b/Stochastisches_Modell.py index bd882bd..4e85267 100644 --- a/Stochastisches_Modell.py +++ b/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) \ No newline at end of file + def update_sigma(self, sigma_hat_dict): + for g, val in sigma_hat_dict.items(): + self.sigma0_groups[g] = float(val) \ No newline at end of file