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 msr::sensitivity::TrackerSensitivity;
25use statrs::distribution::{ContinuousCDF, Normal};
26use std::ops::Add;
27
28use super::ODSolution;
29
30impl<StateType, EstType, MsrSize, Trk> ODSolution<StateType, EstType, MsrSize, Trk>
31where
32    StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType>,
33    EstType: Estimate<StateType>,
34    MsrSize: DimName,
35    Trk: TrackerSensitivity<StateType, StateType>,
36    <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send,
37    DefaultAllocator: Allocator<<StateType as State>::Size>
38        + Allocator<<StateType as State>::VecLength>
39        + Allocator<MsrSize>
40        + Allocator<MsrSize, <StateType as State>::Size>
41        + Allocator<MsrSize, MsrSize>
42        + Allocator<<StateType as State>::Size, <StateType as State>::Size>
43        + Allocator<<StateType as State>::Size, MsrSize>,
44{
45    /// Returns the root mean square of the prefit residuals
46    pub fn rms_prefit_residuals(&self) -> f64 {
47        let mut sum = 0.0;
48        for residual in self.residuals.iter().flatten() {
49            sum += residual.prefit.dot(&residual.prefit);
50        }
51        (sum / (self.residuals.len() as f64)).sqrt()
52    }
53
54    /// Returns the root mean square of the postfit residuals
55    pub fn rms_postfit_residuals(&self) -> f64 {
56        let mut sum = 0.0;
57        for residual in self.residuals.iter().flatten() {
58            sum += residual.postfit.dot(&residual.postfit);
59        }
60        (sum / (self.residuals.len() as f64)).sqrt()
61    }
62
63    /// Returns the root mean square of the prefit residual ratios
64    pub fn rms_residual_ratios(&self) -> f64 {
65        let mut sum = 0.0;
66        for residual in self.residuals.iter().flatten() {
67            sum += residual.ratio.powi(2);
68        }
69        (sum / (self.residuals.len() as f64)).sqrt()
70    }
71
72    /// Computes the fraction of residual ratios that lie within ±threshold.
73    pub fn residual_ratio_within_threshold(&self, threshold: f64) -> Result<f64, ODError> {
74        let mut total = 0;
75        let mut count_within = 0;
76        for residual in self.residuals.iter().flatten() {
77            total += 1;
78            if residual.ratio.abs() <= threshold {
79                count_within += 1;
80            }
81        }
82        if total > 0 {
83            Ok(count_within as f64 / total as f64)
84        } else {
85            Err(ODError::ODNoResiduals {
86                action: "percentage of residuals within threshold",
87            })
88        }
89    }
90
91    /// Computes the Kolmogorov–Smirnov statistic for the aggregated residual ratios, by tracker and measurement type.
92    ///
93    /// Returns Ok(ks_statistic) if residuals are available.
94    pub fn ks_test_normality(&self) -> Result<f64, ODError> {
95        let mut residual_ratios = self
96            .residuals
97            .iter()
98            .flatten()
99            .map(|resid| resid.ratio)
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 residual ratios 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: "percentage of residuals within threshold",
149            });
150        }
151        let ks_stat = self.ks_test_normality()?;
152
153        // Default to 5% probability
154        let alpha = alpha.unwrap_or(0.05);
155
156        // Compute the c_alpha
157        let c_alpha = (-(alpha * 0.5).ln() * 0.5).sqrt();
158
159        let d_critical = c_alpha / (n as f64).sqrt();
160
161        Ok(ks_stat <= d_critical)
162    }
163}
nyx_space/od/process/solution/stats.rs

nyx_space/od/process/solution/
stats.rs
