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, by tracker and measurement type.
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            .residuals
98            .iter()
99            .flatten()
100            .map(|resid| resid.ratio)
101            .collect::<Vec<f64>>();
102
103        if residual_ratios.is_empty() {
104            return Err(ODError::ODNoResiduals {
105                action: "percentage of residuals within threshold",
106            });
107        }
108        residual_ratios.sort_by(|a, b| a.partial_cmp(b).unwrap());
109        let n = residual_ratios.len() as f64;
110        let mean = residual_ratios.iter().sum::<f64>() / n;
111        let variance = residual_ratios
112            .iter()
113            .map(|v| (v - mean).powi(2))
114            .sum::<f64>()
115            / n;
116        let std_dev = variance.sqrt();
117
118        // Create a normal distribution based on the empirical mean and std deviation.
119        let normal_dist = Normal::new(mean, std_dev).unwrap();
120
121        // Compute the maximum difference between the empirical CDF and the normal CDF.
122        let mut d_stat = 0.0;
123        for (i, &value) in residual_ratios.iter().enumerate() {
124            let empirical_cdf = (i + 1) as f64 / n;
125            let model_cdf = normal_dist.cdf(value);
126            let diff = (empirical_cdf - model_cdf).abs();
127            if diff > d_stat {
128                d_stat = diff;
129            }
130        }
131        Ok(d_stat)
132    }
133
134    /// Checks whether the residual ratios pass a normality test at a given significance level `alpha`, default to 0.05.
135    ///
136    /// This uses a simplified KS-test threshold: D_alpha = c(α) / √n.
137    /// For example, for α = 0.05, c(α) is approximately 1.36.
138    /// α=0.05 means a 5% probability of Type I error (incorrectly rejecting the null hypothesis when it is true).
139    /// This threshold is standard in many statistical tests to balance sensitivity and false positives
140    ///
141    /// 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
142    ///
143    /// Returns Ok(true) if the residuals are consistent with a normal distribution,
144    /// Ok(false) if not, or None if no residuals are available.
145    pub fn is_normal(&self, alpha: Option<f64>) -> Result<bool, ODError> {
146        let n = self.residuals.iter().flatten().count();
147        if n == 0 {
148            return Err(ODError::ODNoResiduals {
149                action: "evaluating residual normality",
150            });
151        } else if n < 35 {
152            warn!("KS normality test unreliable for n={n} < 35; skipping");
153        }
154        let ks_stat = self.ks_test_normality()?;
155
156        // Default to 5% probability
157        let alpha = alpha.unwrap_or(0.05);
158
159        // Compute the c_alpha
160        let c_alpha = (-(alpha * 0.5).ln() * 0.5).sqrt();
161
162        let d_critical = c_alpha / (n as f64).sqrt();
163
164        Ok(ks_stat <= d_critical)
165    }
166}
nyx_space/od/process/solution/stats.rs

nyx_space/od/process/solution/
stats.rs
