nyx_space/od/process/solution/
stats.rs

1/*
2    Nyx, blazing fast astrodynamics
3    Copyright (C) 2018-onwards Christopher Rabotin <christopher.rabotin@gmail.com>
4
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU Affero General Public License as published
7    by the Free Software Foundation either version 3 of the License, or
8    (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU Affero General Public License for more details.
14
15    You should have received a copy of the GNU Affero General Public License
16    along with this program.  If not, see <https://www.gnu.org/licenses/>.
17*/
18
19use crate::linalg::allocator::Allocator;
20use crate::linalg::{DefaultAllocator, DimName};
21use crate::md::trajectory::Interpolatable;
22pub use crate::od::estimate::*;
23pub use crate::od::*;
24use log::warn;
25use msr::sensitivity::TrackerSensitivity;
26use statrs::distribution::{ContinuousCDF, Normal};
27use std::ops::Add;
28
29use super::ODSolution;
30
31impl<StateType, EstType, MsrSize, Trk> ODSolution<StateType, EstType, MsrSize, Trk>
32where
33    StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType>,
34    EstType: Estimate<StateType>,
35    MsrSize: DimName,
36    Trk: TrackerSensitivity<StateType, StateType>,
37    <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send,
38    DefaultAllocator: Allocator<<StateType as State>::Size>
39        + Allocator<<StateType as State>::VecLength>
40        + Allocator<MsrSize>
41        + Allocator<MsrSize, <StateType as State>::Size>
42        + Allocator<MsrSize, MsrSize>
43        + Allocator<<StateType as State>::Size, <StateType as State>::Size>
44        + Allocator<<StateType as State>::Size, MsrSize>,
45{
46    /// Returns the root mean square of the prefit residuals
47    pub fn rms_prefit_residuals(&self) -> f64 {
48        let mut sum = 0.0;
49        for residual in self.residuals.iter().flatten() {
50            sum += residual.prefit.dot(&residual.prefit);
51        }
52        (sum / (self.residuals.len() as f64)).sqrt()
53    }
54
55    /// Returns the root mean square of the postfit residuals
56    pub fn rms_postfit_residuals(&self) -> f64 {
57        let mut sum = 0.0;
58        for residual in self.residuals.iter().flatten() {
59            sum += residual.postfit.dot(&residual.postfit);
60        }
61        (sum / (self.residuals.len() as f64)).sqrt()
62    }
63
64    /// Returns the root mean square of the prefit residual ratios
65    pub fn rms_residual_ratios(&self) -> f64 {
66        let mut sum = 0.0;
67        for residual in self.residuals.iter().flatten() {
68            sum += residual.ratio.powi(2);
69        }
70        (sum / (self.residuals.len() as f64)).sqrt()
71    }
72
73    /// Computes the fraction of residual ratios that lie within ±threshold.
74    pub fn residual_ratio_within_threshold(&self, threshold: f64) -> Result<f64, ODError> {
75        let mut total = 0;
76        let mut count_within = 0;
77        for residual in self.residuals.iter().flatten() {
78            total += 1;
79            if residual.ratio.abs() <= threshold {
80                count_within += 1;
81            }
82        }
83        if total > 0 {
84            Ok(count_within as f64 / total as f64)
85        } else {
86            Err(ODError::ODNoResiduals {
87                action: "percentage of residuals within threshold",
88            })
89        }
90    }
91
92    /// Computes the Kolmogorov–Smirnov statistic for the aggregated residual ratios of the accepted residuals.
93    ///
94    /// Returns Ok(ks_statistic) if residuals are available.
95    pub fn ks_test_normality(&self) -> Result<f64, ODError> {
96        let mut residual_ratios = self
97            .accepted_residuals()
98            .iter()
99            .flat_map(|resid| resid.whitened_resid.into_iter().copied())
100            .collect::<Vec<f64>>();
101
102        if residual_ratios.is_empty() {
103            return Err(ODError::ODNoResiduals {
104                action: "percentage of residuals within threshold",
105            });
106        }
107        residual_ratios.sort_by(|a, b| a.partial_cmp(b).unwrap());
108        let n = residual_ratios.len() as f64;
109        let mean = residual_ratios.iter().sum::<f64>() / n;
110        let variance = residual_ratios
111            .iter()
112            .map(|v| (v - mean).powi(2))
113            .sum::<f64>()
114            / n;
115        let std_dev = variance.sqrt();
116
117        // Create a normal distribution based on the empirical mean and std deviation.
118        let normal_dist = Normal::new(mean, std_dev).unwrap();
119
120        // Compute the maximum difference between the empirical CDF and the normal CDF.
121        let mut d_stat = 0.0;
122        for (i, &value) in residual_ratios.iter().enumerate() {
123            let empirical_cdf = (i + 1) as f64 / n;
124            let model_cdf = normal_dist.cdf(value);
125            let diff = (empirical_cdf - model_cdf).abs();
126            if diff > d_stat {
127                d_stat = diff;
128            }
129        }
130        Ok(d_stat)
131    }
132
133    /// Checks whether the whitened residuals of the accepted residuals pass a normality test at a given significance level `alpha`, default to 0.05.
134    ///
135    /// This uses a simplified KS-test threshold: D_alpha = c(α) / √n.
136    /// For example, for α = 0.05, c(α) is approximately 1.36.
137    /// α=0.05 means a 5% probability of Type I error (incorrectly rejecting the null hypothesis when it is true).
138    /// This threshold is standard in many statistical tests to balance sensitivity and false positives
139    ///
140    /// The computation of the c(alpha) is from https://en.wikipedia.org/w/index.php?title=Kolmogorov%E2%80%93Smirnov_test&oldid=1280260701#Two-sample_Kolmogorov%E2%80%93Smirnov_test
141    ///
142    /// Returns Ok(true) if the residuals are consistent with a normal distribution,
143    /// Ok(false) if not, or None if no residuals are available.
144    pub fn is_normal(&self, alpha: Option<f64>) -> Result<bool, ODError> {
145        let n = self.residuals.iter().flatten().count();
146        if n == 0 {
147            return Err(ODError::ODNoResiduals {
148                action: "evaluating residual normality",
149            });
150        } else if n < 35 {
151            warn!("KS normality test unreliable for n={n} < 35; skipping");
152        }
153        let ks_stat = self.ks_test_normality()?;
154
155        // Default to 5% probability
156        let alpha = alpha.unwrap_or(0.05);
157
158        // Compute the c_alpha
159        let c_alpha = (-(alpha * 0.5).ln() * 0.5).sqrt();
160
161        let d_critical = c_alpha / (n as f64).sqrt();
162
163        Ok(ks_stat <= d_critical)
164    }
165
166    /// Checks whether the filter innovations are statistically consistent
167    /// by performing a Chi-squared test on the Normalized Innovation Squared (NIS).
168    ///
169    /// For each accepted residual, NIS is computed as:
170    /// ```text
171    ///     prefit^T * S_k^-1 * prefit
172    /// ```
173    ///
174    /// The sum of NIS values should fall within the confidence interval of a
175    /// Chi-squared distribution with degrees of freedom `k = n * m`, where `n`
176    /// is the number of residuals and `m` is the measurement dimension.
177    ///
178    /// Returns Ok(true) if the filter is consistent, Ok(false) if the filter
179    /// is over-confident or under-confident, or an error if no residuals are available.
180    pub fn is_nis_consistent(&self, alpha: Option<f64>) -> Result<bool, ODError> {
181        let n = self.accepted_residuals().iter().count();
182
183        if n == 0 {
184            return Err(ODError::ODNoResiduals {
185                action: "evaluating NIS consistency",
186            });
187        }
188
189        // Sum the NIS across all residuals.
190        // Assuming your Residual struct has an `nis` field or a method to compute it
191        // from the ratio: `ratio.powi(2) * measurement_dim`
192        let nis_sum: f64 = self.accepted_residuals().iter().map(|r| r.nis()).sum();
193
194        // Total degrees of freedom: number of residuals * measurement dimension.
195        // Adjust `M::DIM` based on how you access the dimension in this scope.
196        let k = (n * MsrSize::DIM) as f64;
197        if k < 35.0 {
198            warn!("NIS consistency test lacks statistical power for n*MsrSize={k} < 35");
199        }
200        // Default to a 5% probability of Type I error (95% confidence interval)
201        let alpha = alpha.unwrap_or(0.05);
202
203        // For a two-sided test, we need the standard normal quantile for 1 - (alpha / 2).
204        // If alpha = 0.05, the critical z-score is approximately 1.95996.
205        let z_critical = Normal::new(0.0, 1.0)
206            .unwrap()
207            .inverse_cdf(1.0 - alpha / 2.0);
208
209        // Use the Wilson-Hilferty transformation to approximate the Chi-squared
210        // lower and upper critical bounds.
211        let factor = 2.0 / (9.0 * k);
212        let lower_bound = k * (1.0 - factor - z_critical * factor.sqrt()).powi(3);
213        let upper_bound = k * (1.0 - factor + z_critical * factor.sqrt()).powi(3);
214
215        if nis_sum > upper_bound {
216            warn!("NIS consistency test failed high: NIS sum {nis_sum:.6} > upper bound {upper_bound:.6}. Innovations are larger than expected.");
217            warn!("Filter may be overconfident: P, R, or Q may be too small, or the dynamics/measurement model may be biased.");
218            Ok(false)
219        } else if nis_sum < lower_bound {
220            warn!("NIS consistency test failed low: NIS sum {nis_sum:.6} < lower bound {lower_bound:.6}. Innovations are smaller than expected.");
221            warn!("Filter may be underconfident: P, R, or Q may be too large.");
222            Ok(false)
223        } else {
224            Ok(true)
225        }
226    }
227}
nyx_space/od/process/solution/stats.rs

nyx_space/od/process/solution/
stats.rs
