From 26059ecee21b6c2f2af005571051b11d669856f7 Mon Sep 17 00:00:00 2001 From: N3uralN3twork Date: Sun, 2 Aug 2026 14:44:32 -0400 Subject: [PATCH 1/2] Solidified numeric policy for crate --- benches/summary_statistics.rs | 24 +- docs/SUMMARY_STATISTICS.md | 197 +++++++++ src/correlation.rs | 186 +++++++- src/deviation.rs | 304 ++++++++++++- src/entropy.rs | 122 +++++- src/errors.rs | 248 +++++++++++ src/lib.rs | 13 +- src/policies.rs | 84 +++- src/summary_statistics/descriptive.rs | 81 ++-- src/summary_statistics/means.rs | 590 +++++++++++++++++++++++++- src/summary_statistics/mod.rs | 235 ++++++++-- tests/numeric_policy.rs | 103 +++++ tests/policies.rs | 32 +- 13 files changed, 2099 insertions(+), 120 deletions(-) create mode 100644 docs/SUMMARY_STATISTICS.md create mode 100644 tests/numeric_policy.rs diff --git a/benches/summary_statistics.rs b/benches/summary_statistics.rs index b1b607b..1f9013a 100644 --- a/benches/summary_statistics.rs +++ b/benches/summary_statistics.rs @@ -4,7 +4,9 @@ use criterion::{ use ndarray::prelude::*; use ndarray_rand::rand_distr::Uniform; use ndarray_rand::RandomExt; -use ndarray_stats::{DescriptiveStatistics, QuantileExt, SummaryStatisticsExt}; +use ndarray_stats::{ + policies::NumericPolicy, DescriptiveStatistics, QuantileExt, SummaryStatisticsExt, +}; mod common; @@ -77,7 +79,10 @@ fn descriptive_statistics(c: &mut Criterion) { group.bench_with_input(format!("fused/{len}"), &data, |b, data| { b.iter(|| { - let summary = black_box(data.descriptive_statistics().unwrap()); + let summary = black_box( + data.descriptive_statistics_with_policy(NumericPolicy::default()) + .unwrap(), + ); black_box(score_f64(&summary)); }) }); @@ -113,7 +118,10 @@ fn descriptive_statistics_axis(c: &mut Criterion) { group.bench_with_input(format!("fused_axis0/{len}"), &data, |b, data| { b.iter(|| { - let summaries = black_box(data.descriptive_statistics_axis(Axis(0)).unwrap()); + let summaries = black_box( + data.descriptive_statistics_axis_with_policy(Axis(0), NumericPolicy::default()) + .unwrap(), + ); let score = summaries.iter().map(score_f64).sum::(); black_box(score); }) @@ -127,7 +135,10 @@ fn descriptive_statistics_axis(c: &mut Criterion) { group.bench_with_input(format!("fused_axis1/{len}"), &data, |b, data| { b.iter(|| { - let summaries = black_box(data.descriptive_statistics_axis(Axis(1)).unwrap()); + let summaries = black_box( + data.descriptive_statistics_axis_with_policy(Axis(1), NumericPolicy::default()) + .unwrap(), + ); let score = summaries.iter().map(score_f64).sum::(); black_box(score); }) @@ -152,7 +163,10 @@ fn descriptive_statistics_f32(c: &mut Criterion) { group.bench_function(format!("fused/{len}"), |b| { let data: Array1 = Array::random(len, Uniform::new(0.0, 1.0).unwrap()); b.iter(|| { - let summary = black_box(data.descriptive_statistics().unwrap()); + let summary = black_box( + data.descriptive_statistics_with_policy(NumericPolicy::default()) + .unwrap(), + ); black_box(score_f32(&summary)); }) }); diff --git a/docs/SUMMARY_STATISTICS.md b/docs/SUMMARY_STATISTICS.md new file mode 100644 index 0000000..8b69789 --- /dev/null +++ b/docs/SUMMARY_STATISTICS.md @@ -0,0 +1,197 @@ +# Summary statistics + +This guide covers the public summary-statistics API implemented in +[`summary_statistics/mod.rs`](src/summary_statistics/mod.rs), +[`summary_statistics/descriptive.rs`](src/summary_statistics/descriptive.rs), +and [`summary_statistics/means.rs`](src/summary_statistics/means.rs). + +The methods are provided by the `SummaryStatisticsExt` trait. Bring the trait +into scope before calling them on an `ndarray` array: + +```rust +use ndarray::{array, Axis}; +use ndarray_stats::{policies::NumericPolicy, SummaryStatisticsExt}; + +let x = array![1.0, 2.0, 3.0, 4.0]; +let weights = array![1.0, 2.0, 1.0, 2.0]; +let matrix = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]; +let axis_weights = array![1.0, 2.0, 1.0]; +let missing = array![1.0, f64::NAN, 3.0]; +let missing_weights = array![1.0, 2.0, 1.0]; +let matrix_with_missing = array![[1.0, f64::NAN, 3.0], [4.0, 5.0, 6.0]]; +``` + +## Important conventions + +- Legacy methods such as `mean` and `descriptive_statistics` return their + historical error types. Policy-aware methods have a `_with_policy` suffix + and return `errors::StatisticalError`. +- `NumericPolicy::propagate()` preserves ordinary floating-point behavior. + `NumericPolicy::omit_missing()` omits `NaN` values, and + `NumericPolicy::reject_non_finite()` returns an error for `NaN` or infinity. +- Axis methods reduce each lane along the selected axis. The result has the + input shape with that axis removed. +- For moment vectors, element `i` is the moment of order `i`; therefore + `raw_moments(3)` returns orders `0` through `3`. +- For weighted variance and standard deviation, `ddof = 0.0` is the + population calculation and `ddof = 1.0` is the sample calculation. + +## `DescriptiveStatistics` + +`DescriptiveStatistics` is a compact summary containing the count, mean, +minimum, maximum, and variance accumulator for a non-empty input. Create one +with `descriptive_statistics` or an axis variant. + +```rust +let summary = x.descriptive_statistics().unwrap(); +``` + +| Method | Purpose | Example result | +| --- | --- | ---: | +| `count()` | Number of included observations. | `summary.count()` -> `4` | +| `mean()` | Arithmetic mean. | `summary.mean()` -> `2.5` | +| `min()` | Smallest observation. | `summary.min()` -> `1.0` | +| `max()` | Largest observation. | `summary.max()` -> `4.0` | +| `population_variance()` | Variance divided by `n`. | `summary.population_variance()` -> `1.25` | +| `sample_variance()` | Variance divided by `n - 1`. | `summary.sample_variance()` -> `1.666...` | +| `population_std()` | Square root of population variance. | `summary.population_std()` -> `1.118...` | +| `sample_std()` | Square root of sample variance. | `summary.sample_std()` -> `1.291...` | + +For a one-observation summary, the sample variance and sample standard +deviation follow floating-point division semantics and are `NaN`. + +## Descriptive summary methods + +| Method | Purpose | Example | +| --- | --- | --- | +| `descriptive_statistics()` | Summarize all elements using the compatibility policy. | `x.descriptive_statistics().unwrap().mean()` -> `2.5` | +| `descriptive_statistics_with_policy(policy)` | Summarize all elements with an explicit numeric policy. | `missing.descriptive_statistics_with_policy(NumericPolicy::omit_missing()).unwrap().count()` -> `2` | +| `descriptive_statistics_axis(axis)` | Return one summary for every lane along an axis. | `matrix.descriptive_statistics_axis(Axis(1)).unwrap()[0].mean()` -> `2.0` | +| `descriptive_statistics_axis_with_policy(axis, policy)` | Return per-lane summaries with an explicit policy. | `matrix_with_missing.descriptive_statistics_axis_with_policy(Axis(1), NumericPolicy::omit_missing()).unwrap()[0].count()` -> `2` | + +## Means and weighted reductions + +| Method | Purpose | Example result | +| --- | --- | ---: | +| `mean()` | Arithmetic mean of all elements. | `x.mean().unwrap()` -> `2.5` | +| `weighted_mean(weights)` | Mean using the supplied weights. | `x.weighted_mean(&weights).unwrap()` -> `2.666...` | +| `weighted_sum(weights)` | Sum of `value * weight` pairs. | `x.weighted_sum(&weights).unwrap()` -> `16.0` | +| `weighted_mean_axis(axis, weights)` | Weighted mean of every lane. | `matrix.weighted_mean_axis(Axis(1), &axis_weights).unwrap()` -> `[2.0, 5.0]` | +| `weighted_sum_axis(axis, weights)` | Weighted sum of every lane. | `matrix.weighted_sum_axis(Axis(1), &axis_weights).unwrap()` -> `[8.0, 20.0]` | +| `harmonic_mean()` | Harmonic mean. | `x.harmonic_mean().unwrap()` -> `1.92` | +| `geometric_mean()` | Geometric mean. | `x.geometric_mean().unwrap()` -> `24^(1/4) ~ 2.213` | + +Weighted reductions require matching shapes. Axis-weight arrays must have one +element for every value in the selected axis. + +## Modes + +Modes use `PartialEq`, so ordinary floating-point arrays can use these methods +without `Hash` or `Eq`. Ties preserve first-occurrence order. + +```rust +let mode_input = array![1, 2, 2, 3, 3]; +let mode_matrix = array![[1, 2, 2], [3, 3, 4]]; +``` + +| Method | Purpose | Example result | +| --- | --- | --- | +| `mode()` | Return the first value among the modes. | `mode_input.mode().unwrap()` -> `2` | +| `modes()` | Return all modes in first-occurrence order. | `mode_input.modes().unwrap()` -> `[2, 3]` | +| `mode_axis(axis)` | Return the first mode for every lane. | `mode_matrix.mode_axis(Axis(1)).unwrap()` -> `[2, 3]` | + +## Moments and shape statistics + +| Method | Purpose | Example result | +| --- | --- | ---: | +| `raw_moment(order)` | Return one raw moment, `mean(x^order)`. | `x.raw_moment(2).unwrap()` -> `7.5` | +| `raw_moments(order)` | Return raw moments from order `0` through `order`. | `x.raw_moments(3).unwrap()` -> `[1.0, 2.5, 7.5, 25.0]` | +| `central_moment(order)` | Return one central moment, `mean((x - mean)^order)`. | `x.central_moment(2).unwrap()` -> `1.25` | +| `central_moments(order)` | Return central moments from order `0` through `order`. | `x.central_moments(3).unwrap()` -> `[1.0, 0.0, 1.25, 0.0]` | +| `standardized_moment(order)` | Return one central moment divided by the corresponding power of standard deviation. | `x.standardized_moment(3).unwrap()` -> `0.0` | +| `standardized_moments(order)` | Return standardized moments from order `0` through `order`. | `x.standardized_moments(4).unwrap()` -> `[1.0, 0.0, 1.0, 0.0, 1.64]` | +| `skewness()` | Return the third standardized moment. | `x.skewness().unwrap()` -> `0.0` | +| `kurtosis()` | Return Pearson's kurtosis, the fourth standardized moment. | `x.kurtosis().unwrap()` -> `1.64` | + +The zeroth raw and standardized moments are `1.0`; the first central and +standardized moments are `0.0`. A zero-variance input generally produces `NaN` +for standardized moments of order two or greater. + +## Weighted variance and standard deviation + +| Method | Purpose | Example result | +| --- | --- | ---: | +| `weighted_var(weights, ddof)` | Weighted variance for all elements. | `x.weighted_var(&weights, 0.0).unwrap()` -> `1.222...` | +| `weighted_std(weights, ddof)` | Square root of weighted variance. | `x.weighted_std(&weights, 0.0).unwrap()` -> `1.105...` | +| `weighted_var_axis(axis, weights, ddof)` | Weighted variance for every lane. | `matrix.weighted_var_axis(Axis(1), &axis_weights, 0.0).unwrap()` -> `[0.5, 0.5]` | +| `weighted_std_axis(axis, weights, ddof)` | Weighted standard deviation for every lane. | `matrix.weighted_std_axis(Axis(1), &axis_weights, 0.0).unwrap()` -> `[0.707..., 0.707...]` | + +`ddof` must be between `0.0` and `1.0`. A value outside that interval is a +programming error for the legacy methods and may panic. + +## Policy-aware methods + +Policy-aware methods use the same calculations as the legacy methods but make +missing-value and infinity behavior explicit. The examples below use +`NumericPolicy::omit_missing()`, so the `NaN` value is omitted. For paired +operations, the value and its corresponding weight are omitted together. + +### Policy-aware means and reductions + +| Method | Example result | +| --- | ---: | +| `mean_with_policy(policy)` | `missing.mean_with_policy(NumericPolicy::omit_missing()).unwrap()` -> `2.0` | +| `weighted_mean_with_policy(weights, policy)` | `missing.weighted_mean_with_policy(&missing_weights, NumericPolicy::omit_missing()).unwrap()` -> `2.0` | +| `weighted_sum_with_policy(weights, policy)` | `missing.weighted_sum_with_policy(&missing_weights, NumericPolicy::omit_missing()).unwrap()` -> `4.0` | +| `weighted_mean_axis_with_policy(axis, weights, policy)` | `matrix_with_missing.weighted_mean_axis_with_policy(Axis(1), &axis_weights, NumericPolicy::omit_missing()).unwrap()` -> `[2.0, 5.0]` | +| `weighted_sum_axis_with_policy(axis, weights, policy)` | `matrix_with_missing.weighted_sum_axis_with_policy(Axis(1), &axis_weights, NumericPolicy::omit_missing()).unwrap()` -> `[4.0, 20.0]` | +| `harmonic_mean_with_policy(policy)` | Intended harmonic mean after omission: `missing.harmonic_mean_with_policy(NumericPolicy::omit_missing()).unwrap()` -> `1.5` | +| `geometric_mean_with_policy(policy)` | `missing.geometric_mean_with_policy(NumericPolicy::omit_missing()).unwrap()` -> `sqrt(3) ~ 1.732` | + +### Policy-aware moments and shape statistics + +| Method | Example result | +| --- | ---: | +| `raw_moment_with_policy(order, policy)` | `missing.raw_moment_with_policy(2, NumericPolicy::omit_missing()).unwrap()` -> `5.0` | +| `raw_moments_with_policy(order, policy)` | `missing.raw_moments_with_policy(2, NumericPolicy::omit_missing()).unwrap()` -> `[1.0, 2.0, 5.0]` | +| `central_moment_with_policy(order, policy)` | `missing.central_moment_with_policy(2, NumericPolicy::omit_missing()).unwrap()` -> `1.0` | +| `central_moments_with_policy(order, policy)` | `missing.central_moments_with_policy(3, NumericPolicy::omit_missing()).unwrap()` -> `[1.0, 0.0, 1.0, 0.0]` | +| `standardized_moment_with_policy(order, policy)` | `missing.standardized_moment_with_policy(3, NumericPolicy::omit_missing()).unwrap()` -> `0.0` | +| `standardized_moments_with_policy(order, policy)` | `missing.standardized_moments_with_policy(3, NumericPolicy::omit_missing()).unwrap()` -> `[1.0, 0.0, 1.0, 0.0]` | +| `skewness_with_policy(policy)` | `missing.skewness_with_policy(NumericPolicy::omit_missing()).unwrap()` -> `0.0` | +| `kurtosis_with_policy(policy)` | `missing.kurtosis_with_policy(NumericPolicy::omit_missing()).unwrap()` -> `1.0` | + +### Policy-aware weighted variance + +| Method | Example result | +| --- | ---: | +| `weighted_var_with_policy(weights, ddof, policy)` | `missing.weighted_var_with_policy(&missing_weights, 0.0, NumericPolicy::omit_missing()).unwrap()` -> `1.0` | +| `weighted_std_with_policy(weights, ddof, policy)` | `missing.weighted_std_with_policy(&missing_weights, 0.0, NumericPolicy::omit_missing()).unwrap()` -> `1.0` | +| `weighted_var_axis_with_policy(axis, weights, ddof, policy)` | `matrix_with_missing.weighted_var_axis_with_policy(Axis(1), &axis_weights, 0.0, NumericPolicy::omit_missing()).unwrap()` -> `[1.0, 0.5]` | +| `weighted_std_axis_with_policy(axis, weights, ddof, policy)` | `matrix_with_missing.weighted_std_axis_with_policy(Axis(1), &axis_weights, 0.0, NumericPolicy::omit_missing()).unwrap()` -> `[1.0, 0.707...]` | + +With `NumericPolicy::reject_non_finite()`, the same calls return a +`StatisticalError` when they encounter `NaN` or infinity. With +`NumericPolicy::propagate()`, no filtering occurs and IEEE-754 behavior is +preserved. + +The current `harmonic_mean_with_policy` implementation returns the reciprocal +of the intended harmonic mean (`0.666...` for `missing`), while the legacy +`harmonic_mean` returns `1.5`. The table documents the method's stated +statistical meaning; this implementation discrepancy is not changed by this +documentation-only update. + +## Errors and edge cases + +- Empty inputs return `EmptyInput` for legacy scalar methods and an + `EmptyInput` variant of `StatisticalError` for policy-aware methods. +- Weighted operations return a shape-mismatch error when values and weights do + not align. +- Axis methods panic when the requested axis is out of bounds. +- Moment orders are `u16`; very large orders can overflow the internal `i32` + power representation. +- `sample_variance()` and `sample_std()` are undefined for a single + observation and return `NaN` according to floating-point semantics. + +Private accumulators and helper functions in the implementation modules are +not part of the public API and are therefore not listed here. diff --git a/src/correlation.rs b/src/correlation.rs index 93d0e50..f89d22b 100644 --- a/src/correlation.rs +++ b/src/correlation.rs @@ -1,4 +1,5 @@ -use crate::errors::EmptyInput; +use crate::errors::{EmptyInput, StatisticalError, StatisticalErrorContext, StatisticalOperation}; +use crate::policies::{apply_numeric_policy, NumericPolicy}; use ndarray::prelude::*; use num_traits::{Float, FromPrimitive}; use std::cmp::Ordering; @@ -159,6 +160,36 @@ pub trait CorrelationExt { where A: Float + FromPrimitive; + /// Policy-aware covariance. Missing observations are omitted listwise: + /// an observation is retained only when every variable has a usable value. + fn cov_with_policy( + &self, + ddof: A, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive; + + /// Policy-aware Pearson correlation using the rows-as-variables, + /// columns-as-observations convention of [`cov_with_policy`](Self::cov_with_policy). + fn pearson_correlation_with_policy( + &self, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive; + + fn spearman_correlation_with_policy( + &self, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive; + + fn kendall_tau_with_policy(&self, policy: NumericPolicy) -> Result, StatisticalError> + where + A: Float + FromPrimitive; + private_decl! {} } @@ -241,9 +272,162 @@ impl CorrelationExt for ArrayRef2 { Ok(result) } + fn cov_with_policy(&self, ddof: A, policy: NumericPolicy) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + { + covariance_with_policy(self, ddof, policy, StatisticalOperation::Covariance) + } + + fn pearson_correlation_with_policy( + &self, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::PearsonCorrelation); + let (n_variables, n_observations) = self.dim(); + if n_variables == 0 || n_observations == 0 { + return Err(StatisticalError::EmptyInput { context }); + } + let covariance = covariance_with_policy( + self, + A::zero(), + policy, + StatisticalOperation::PearsonCorrelation, + )?; + let standard_deviations = + Array1::from_iter((0..n_variables).map(|index| covariance[[index, index]].sqrt())); + let std_matrix = standard_deviations + .view() + .insert_axis(Axis(1)) + .dot(&standard_deviations.view().insert_axis(Axis(0))); + Ok(covariance / std_matrix) + } + + fn spearman_correlation_with_policy( + &self, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + { + let operation = StatisticalOperation::SpearmanCorrelation; + let ranks = policy_matrix(self, policy, operation)?; + let ranks = rank_rows(&ranks).map_err(|_| StatisticalError::EmptyInput { + context: StatisticalErrorContext::new(operation), + })?; + ranks.pearson_correlation_with_policy(NumericPolicy::propagate()) + } + + fn kendall_tau_with_policy(&self, policy: NumericPolicy) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + { + let operation = StatisticalOperation::KendallTau; + let data = policy_matrix(self, policy, operation)?; + let (n_variables, n_observations) = data.dim(); + if n_variables == 0 || n_observations == 0 { + return Err(StatisticalError::EmptyInput { + context: StatisticalErrorContext::new(operation), + }); + } + + let mut result = Array2::from_elem((n_variables, n_variables), A::nan()); + for first in 0..n_variables { + for second in first..n_variables { + let tau = kendall_tau_pair(data.row(first), data.row(second)); + result[[first, second]] = tau; + result[[second, first]] = tau; + } + } + Ok(result) + } + private_impl! {} } +fn policy_matrix( + data: &ArrayRef2, + policy: NumericPolicy, + operation: StatisticalOperation, +) -> Result, StatisticalError> +where + A: Float, +{ + let (n_variables, n_observations) = data.dim(); + let context = StatisticalErrorContext::new(operation); + if n_observations == 0 { + return Err(StatisticalError::EmptyInput { context }); + } + if n_variables == 0 { + return Ok(Array2::zeros((0, n_observations))); + } + + let mut filtered = vec![Vec::new(); n_variables]; + for observation in 0..n_observations { + let mut values = Vec::with_capacity(n_variables); + let mut include = true; + for variable in 0..n_variables { + let context = + StatisticalErrorContext::axis_lane(operation, 1, variable).with_index(observation); + match apply_numeric_policy(data[[variable, observation]], policy, context)? { + Some(value) => values.push(value), + None => { + include = false; + values.push(A::zero()); + } + } + } + if include { + for variable in 0..n_variables { + filtered[variable].push(values[variable]); + } + } + } + + let retained = filtered[0].len(); + if retained == 0 { + return Err(StatisticalError::EmptyInput { context }); + } + let values = filtered + .into_iter() + .flat_map(|row| row.into_iter()) + .collect(); + Ok(Array2::from_shape_vec((n_variables, retained), values) + .expect("filtered covariance rows must have equal lengths")) +} + +fn covariance_with_policy( + data: &ArrayRef2, + ddof: A, + policy: NumericPolicy, + operation: StatisticalOperation, +) -> Result, StatisticalError> +where + A: Float + FromPrimitive + 'static, +{ + let context = StatisticalErrorContext::new(operation); + let filtered = policy_matrix(data, policy, operation)?; + let (n_variables, n_observations) = filtered.dim(); + let n_observations_as_a = A::from_usize(n_observations).expect("length must convert to A"); + if ddof < A::zero() || ddof >= n_observations_as_a { + return Err(StatisticalError::InvalidDegreesOfFreedom { context }); + } + if n_variables == 0 { + return Ok(Array2::zeros((0, 0))); + } + let mean = filtered + .mean_axis(Axis(1)) + .ok_or(StatisticalError::EmptyInput { context })?; + let denoised = &filtered - &mean.insert_axis(Axis(1)); + let denominator = n_observations_as_a - ddof; + Ok(denoised + .dot(&denoised.t()) + .mapv_into(|value| value / denominator)) +} + /// Compute average one-based ranks for every row without modifying the input. fn rank_rows(data: &ArrayRef2) -> Result, EmptyInput> where diff --git a/src/deviation.rs b/src/deviation.rs index 3c35746..5c30a2a 100644 --- a/src/deviation.rs +++ b/src/deviation.rs @@ -1,9 +1,12 @@ use ndarray::{ArrayRef, Dimension, Zip}; -use num_traits::{Signed, ToPrimitive}; +use num_traits::{Float, Signed, ToPrimitive}; use std::convert::Into; use std::ops::AddAssign; -use crate::errors::MultiInputError; +use crate::errors::{ + MultiInputError, StatisticalError, StatisticalErrorContext, StatisticalOperation, +}; +use crate::policies::{apply_numeric_policy_pair, NumericPolicy}; /// An extension trait for `ndarray` providing functions /// to compute different deviation measures. @@ -203,6 +206,89 @@ where where A: AddAssign + Clone + Signed + ToPrimitive; + /// Policy-aware deviation operations. With omission, a pair is omitted + /// when either input value is NaN. + fn count_eq_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + + fn count_neq_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + + fn sq_l2_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + + fn l2_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + + fn l1_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + + fn linf_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + + fn mean_abs_err_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + + fn mean_sq_err_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + + fn root_mean_sq_err_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + + fn peak_signal_to_noise_ratio_with_policy( + &self, + other: &ArrayRef, + maxv: A, + policy: NumericPolicy, + ) -> Result + where + A: Float; + private_decl! {} } @@ -351,5 +437,219 @@ where Ok(psnr) } + fn count_eq_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let pairs = collect_policy_pairs( + self, + other, + policy, + StatisticalErrorContext::new(StatisticalOperation::Deviation), + )?; + Ok(pairs.iter().filter(|(a, b)| a == b).count()) + } + + fn count_neq_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let pairs = collect_policy_pairs( + self, + other, + policy, + StatisticalErrorContext::new(StatisticalOperation::Deviation), + )?; + Ok(pairs.iter().filter(|(a, b)| a != b).count()) + } + + fn sq_l2_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); + let pairs = collect_policy_pairs(self, other, policy, context)?; + if pairs.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(pairs.iter().copied().fold(A::zero(), |sum, (a, b)| { + let difference = a - b; + sum + difference * difference + })) + } + + fn l2_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + Ok(self + .sq_l2_dist_with_policy(other, policy)? + .to_f64() + .expect("failed cast from type A to f64") + .sqrt()) + } + + fn l1_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); + let pairs = collect_policy_pairs(self, other, policy, context)?; + if pairs.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(pairs + .iter() + .copied() + .fold(A::zero(), |sum, (a, b)| sum + (a - b).abs())) + } + + fn linf_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); + let pairs = collect_policy_pairs(self, other, policy, context)?; + if pairs.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(pairs + .iter() + .copied() + .map(|(a, b)| (a - b).abs()) + .fold( + A::zero(), + |max, value| if value > max { value } else { max }, + )) + } + + fn mean_abs_err_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); + let pairs = collect_policy_pairs(self, other, policy, context)?; + if pairs.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + let sum = pairs + .iter() + .copied() + .fold(A::zero(), |sum, (a, b)| sum + (a - b).abs()) + .to_f64() + .expect("failed cast from type A to f64"); + Ok(sum / pairs.len() as f64) + } + + fn mean_sq_err_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let sum = self + .sq_l2_dist_with_policy(other, policy)? + .to_f64() + .expect("failed cast from type A to f64"); + let pairs = collect_policy_pairs( + self, + other, + policy, + StatisticalErrorContext::new(StatisticalOperation::Deviation), + )?; + Ok(sum / pairs.len() as f64) + } + + fn root_mean_sq_err_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + Ok(self.mean_sq_err_with_policy(other, policy)?.sqrt()) + } + + fn peak_signal_to_noise_ratio_with_policy( + &self, + other: &ArrayRef, + maxv: A, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let maxv = maxv.to_f64().expect("failed cast from type A to f64"); + let mse = self.mean_sq_err_with_policy(other, policy)?; + Ok(10. * f64::log10(maxv * maxv / mse)) + } + private_impl! {} } + +fn collect_policy_pairs( + first: &ArrayRef, + second: &ArrayRef, + policy: NumericPolicy, + context: StatisticalErrorContext, +) -> Result, StatisticalError> +where + A: Float, + D: Dimension, +{ + if first.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + if first.shape() != second.shape() { + return Err(StatisticalError::ShapeMismatch { + context, + first_shape: first.shape().to_vec(), + second_shape: second.shape().to_vec(), + }); + } + + first + .iter() + .copied() + .zip(second.iter().copied()) + .enumerate() + .try_fold(Vec::new(), |mut pairs, (index, (first, second))| { + if let Some(pair) = + apply_numeric_policy_pair(first, second, policy, context.with_index(index))? + { + pairs.push(pair); + } + Ok(pairs) + }) +} diff --git a/src/entropy.rs b/src/entropy.rs index 4ba9972..b12b967 100644 --- a/src/entropy.rs +++ b/src/entropy.rs @@ -1,5 +1,9 @@ //! Information theory (e.g. entropy, KL divergence, etc.). -use crate::errors::{EmptyInput, MultiInputError, ShapeMismatch}; +use crate::errors::{ + EmptyInput, MultiInputError, ShapeMismatch, StatisticalError, StatisticalErrorContext, + StatisticalOperation, +}; +use crate::policies::{apply_numeric_policy_pair, collect_with_numeric_policy, NumericPolicy}; use ndarray::{Array, ArrayRef, Dimension, Zip}; use num_traits::Float; @@ -118,6 +122,29 @@ where where A: Float; + /// Policy-aware entropy using the shared missing/non-finite contract. + fn entropy_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float; + + /// Policy-aware KL divergence. Missing pairs are omitted together. + fn kl_divergence_with_policy( + &self, + q: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + + /// Policy-aware cross entropy. Missing pairs are omitted together. + fn cross_entropy_with_policy( + &self, + q: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float; + private_decl! {} } @@ -209,9 +236,102 @@ where Ok(cross_entropy) } + fn entropy_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::Entropy); + let values = collect_with_numeric_policy(self.iter().copied(), policy, context)?; + if values.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(-values.iter().copied().fold(A::zero(), |sum, value| { + if value == A::zero() { + sum + } else { + sum + value * value.ln() + } + })) + } + + fn kl_divergence_with_policy( + &self, + q: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::KLDivergence); + let pairs = collect_policy_pairs(self, q, policy, context)?; + if pairs.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(-pairs.iter().copied().fold(A::zero(), |sum, (p, q)| { + if p == A::zero() { + sum + } else { + sum + p * (q / p).ln() + } + })) + } + + fn cross_entropy_with_policy( + &self, + q: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::CrossEntropy); + let pairs = collect_policy_pairs(self, q, policy, context)?; + if pairs.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(-pairs.iter().copied().fold(A::zero(), |sum, (p, q)| { + if p == A::zero() { + sum + } else { + sum + p * q.ln() + } + })) + } + private_impl! {} } +fn collect_policy_pairs( + p: &ArrayRef, + q: &ArrayRef, + policy: NumericPolicy, + context: StatisticalErrorContext, +) -> Result, StatisticalError> +where + A: Float, + D: Dimension, +{ + if p.shape() != q.shape() { + return Err(StatisticalError::ShapeMismatch { + context, + first_shape: p.shape().to_vec(), + second_shape: q.shape().to_vec(), + }); + } + + p.iter() + .copied() + .zip(q.iter().copied()) + .enumerate() + .try_fold(Vec::new(), |mut pairs, (index, (p, q))| { + if let Some(pair) = apply_numeric_policy_pair(p, q, policy, context.with_index(index))? + { + pairs.push(pair); + } + Ok(pairs) + }) +} + #[cfg(test)] mod tests { use super::EntropyExt; diff --git a/src/errors.rs b/src/errors.rs index 350bd1f..005a34a 100644 --- a/src/errors.rs +++ b/src/errors.rs @@ -27,6 +27,254 @@ pub enum NonFiniteValue { NegativeInfinity, } +/// Identifies the statistical operation that produced an error. +/// +/// The legacy extension methods continue to return their historical error +/// types. Policy-aware methods use [`StatisticalError`] so callers can retain +/// the operation and location that caused a failure. +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub enum StatisticalOperation { + /// The operation is not known, usually because a legacy error was lifted + /// into the shared error model. + Unknown, + DescriptiveStatistics, + Mean, + WeightedMean, + WeightedSum, + WeightedVariance, + WeightedStandardDeviation, + HarmonicMean, + GeometricMean, + RawMoment, + CentralMoment, + StandardizedMoment, + Covariance, + PearsonCorrelation, + SpearmanCorrelation, + KendallTau, + Entropy, + KLDivergence, + CrossEntropy, + Deviation, +} + +impl fmt::Display for StatisticalOperation { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + let name = match self { + StatisticalOperation::Unknown => "statistical operation", + StatisticalOperation::DescriptiveStatistics => "descriptive statistics", + StatisticalOperation::Mean => "mean", + StatisticalOperation::WeightedMean => "weighted mean", + StatisticalOperation::WeightedSum => "weighted sum", + StatisticalOperation::WeightedVariance => "weighted variance", + StatisticalOperation::WeightedStandardDeviation => "weighted standard deviation", + StatisticalOperation::HarmonicMean => "harmonic mean", + StatisticalOperation::GeometricMean => "geometric mean", + StatisticalOperation::RawMoment => "raw moment", + StatisticalOperation::CentralMoment => "central moment", + StatisticalOperation::StandardizedMoment => "standardized moment", + StatisticalOperation::Covariance => "covariance", + StatisticalOperation::PearsonCorrelation => "Pearson correlation", + StatisticalOperation::SpearmanCorrelation => "Spearman correlation", + StatisticalOperation::KendallTau => "Kendall tau", + StatisticalOperation::Entropy => "entropy", + StatisticalOperation::KLDivergence => "KL divergence", + StatisticalOperation::CrossEntropy => "cross entropy", + StatisticalOperation::Deviation => "deviation", + }; + write!(f, "{name}") + } +} + +/// Location metadata attached to a [`StatisticalError`]. +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub struct StatisticalErrorContext { + /// The operation being evaluated. + pub operation: StatisticalOperation, + /// The ndarray axis, when the operation is axis-aware. + pub axis: Option, + /// The zero-based output lane, when the operation is axis-aware. + pub lane: Option, + /// The original zero-based input index within the operation or lane. + pub index: Option, +} + +impl StatisticalErrorContext { + /// Creates context for a whole-array operation. + pub const fn new(operation: StatisticalOperation) -> Self { + Self { + operation, + axis: None, + lane: None, + index: None, + } + } + + /// Creates context for an axis/lane operation. + pub const fn axis_lane(operation: StatisticalOperation, axis: usize, lane: usize) -> Self { + Self { + operation, + axis: Some(axis), + lane: Some(lane), + index: None, + } + } + + /// Creates context for an axis-aware operation before a lane exists. + pub const fn axis(operation: StatisticalOperation, axis: usize) -> Self { + Self { + operation, + axis: Some(axis), + lane: None, + index: None, + } + } + + /// Returns this context with an original input index attached. + pub const fn with_index(self, index: usize) -> Self { + Self { + index: Some(index), + ..self + } + } +} + +impl fmt::Display for StatisticalErrorContext { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "{}", self.operation)?; + if let Some(axis) = self.axis { + write!(f, ", axis {axis}")?; + } + if let Some(lane) = self.lane { + write!(f, ", lane {lane}")?; + } + if let Some(index) = self.index { + write!(f, ", index {index}")?; + } + Ok(()) + } +} + +/// Shared error model for policy-aware statistical operations. +#[derive(Clone, Debug, PartialEq)] +pub enum StatisticalError { + /// No usable observations remained after applying the numeric policy. + EmptyInput { + /// Operation and location metadata. + context: StatisticalErrorContext, + }, + /// A numeric policy rejected a NaN or infinity. + NonFiniteValue { + /// Operation and location metadata. + context: StatisticalErrorContext, + /// The rejected value kind. + value: NonFiniteValue, + }, + /// Two inputs had incompatible shapes. + ShapeMismatch { + /// Operation metadata. + context: StatisticalErrorContext, + /// Shape of the first input. + first_shape: Vec, + /// Shape of the second input. + second_shape: Vec, + }, + /// A required ordering comparison was undefined. + UndefinedOrder { + /// Operation and location metadata. + context: StatisticalErrorContext, + }, + /// A degrees-of-freedom value cannot be used for the available data. + InvalidDegreesOfFreedom { + /// Operation and location metadata. + context: StatisticalErrorContext, + }, + /// The input is outside the operation's mathematical domain. + InvalidDomain { + /// Operation and location metadata. + context: StatisticalErrorContext, + }, +} + +impl StatisticalError { + /// Returns the operation and location attached to this error. + pub fn context(&self) -> StatisticalErrorContext { + match self { + StatisticalError::EmptyInput { context } + | StatisticalError::NonFiniteValue { context, .. } + | StatisticalError::ShapeMismatch { context, .. } + | StatisticalError::UndefinedOrder { context } + | StatisticalError::InvalidDegreesOfFreedom { context } + | StatisticalError::InvalidDomain { context } => *context, + } + } + + /// Returns the operation that produced this error. + pub fn operation(&self) -> StatisticalOperation { + self.context().operation + } +} + +impl fmt::Display for StatisticalError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + let context = self.context(); + match self { + StatisticalError::EmptyInput { .. } => write!(f, "{context}: empty input"), + StatisticalError::NonFiniteValue { value, .. } => write!( + f, + "{context}: {value} is not permitted by the numeric policy" + ), + StatisticalError::ShapeMismatch { + first_shape, + second_shape, + .. + } => write!( + f, + "{context}: array shapes do not match: {first_shape:?} and {second_shape:?}" + ), + StatisticalError::UndefinedOrder { .. } => write!( + f, + "{context}: undefined ordering between a tested pair of values" + ), + StatisticalError::InvalidDegreesOfFreedom { .. } => { + write!(f, "{context}: invalid degrees of freedom") + } + StatisticalError::InvalidDomain { .. } => { + write!(f, "{context}: input is outside the operation's domain") + } + } + } +} + +impl Error for StatisticalError {} + +impl From for StatisticalError { + fn from(_: EmptyInput) -> Self { + StatisticalError::EmptyInput { + context: StatisticalErrorContext::new(StatisticalOperation::Unknown), + } + } +} + +impl From for StatisticalError { + fn from(error: ShapeMismatch) -> Self { + StatisticalError::ShapeMismatch { + context: StatisticalErrorContext::new(StatisticalOperation::Unknown), + first_shape: error.first_shape, + second_shape: error.second_shape, + } + } +} + +impl From for StatisticalError { + fn from(error: MultiInputError) -> Self { + match error { + MultiInputError::EmptyInput => EmptyInput.into(), + MultiInputError::ShapeMismatch(error) => error.into(), + } + } +} + impl fmt::Display for NonFiniteValue { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { match self { diff --git a/src/lib.rs b/src/lib.rs index c935d7c..4acad33 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -10,10 +10,11 @@ //! - [measures of deviation] (count equal, L1, L2 distances, mean squared err etc.) //! - [histogram computation]. //! -//! Numeric policy-aware summary methods are available through -//! [`MissingDataPolicy`], [`InfinityPolicy`], and [`NumericPolicy`]. Existing -//! methods preserve their compatibility behavior, while omission can be -//! requested explicitly with `NumericPolicy::omit_missing()`. +//! Numeric policy-aware methods are available through +//! [`policies::MissingDataPolicy`], [`policies::InfinityPolicy`], and the +//! canonical [`policies::NumericPolicy`] path. Existing methods preserve +//! their compatibility behavior; policy-aware counterparts use a +//! `_with_policy` suffix and the shared [`errors::StatisticalError`] model. //! //! Please feel free to contribute new functionality! A roadmap can be found [here]. //! @@ -39,7 +40,9 @@ pub use crate::deviation::DeviationExt; pub use crate::entropy::EntropyExt; pub use crate::histogram::HistogramExt; pub use crate::maybe_nan::{MaybeNan, MaybeNanExt}; -pub use crate::policies::{InfinityPolicy, MissingDataPolicy, NumericPolicy}; +#[doc(hidden)] +pub use crate::policies::NumericPolicy; +pub use crate::policies::{InfinityPolicy, MissingDataPolicy}; pub use crate::quantile::{interpolate, Quantile1dExt, QuantileExt}; pub use crate::sort::Sort1dExt; pub use crate::summary_statistics::{DescriptiveStatistics, SummaryStatisticsExt}; diff --git a/src/policies.rs b/src/policies.rs index 686a8fa..247a449 100644 --- a/src/policies.rs +++ b/src/policies.rs @@ -1,7 +1,8 @@ //! Shared policies for missing and non-finite numeric values. //! -//! The initial policy surface is used by the policy-aware descriptive-summary -//! methods. It deliberately keeps missingness and infinity separate: +//! The policy surface is shared by policy-aware summary, covariance, +//! correlation, entropy, and deviation methods. It deliberately keeps +//! missingness and infinity separate: //! //! | Input | Policy | Result | //! | --- | --- | --- | @@ -15,9 +16,10 @@ //! No policy imputes, clips, or silently converts a value. In particular, //! omission means omission of `NaN` values only; use [`InfinityPolicy::Reject`] //! when a finite-only calculation is required. Boolean mask workflows and -//! pairwise/listwise deletion for multivariate operations remain follow-up -//! API work, while the existing `*_skipnan` methods remain available for -//! compatibility. +//! Policy-aware two-input operations omit pairs together, and covariance and +//! correlation use listwise omission across variables. Mask-based workflows +//! remain follow-up API work, while the existing `*_skipnan` methods remain +//! available for compatibility. //! //! A finite input is accepted, but finite inputs do not guarantee a finite //! result: arithmetic can overflow or become indeterminate, and those results @@ -25,6 +27,9 @@ //! value filtering occurs; for an order-dependent result, a NaN can therefore //! produce the existing `UndefinedOrder` error rather than a scalar result. +use crate::errors::{NonFiniteValue, StatisticalError, StatisticalErrorContext}; +use num_traits::Float; + /// Controls how `NaN` values are handled when they represent missing data. #[derive(Clone, Copy, Debug, Eq, PartialEq)] pub enum MissingDataPolicy { @@ -102,3 +107,72 @@ impl Default for NumericPolicy { Self::propagate() } } + +/// Applies the shared numeric policy to one value. +pub(crate) fn apply_numeric_policy( + value: A, + policy: NumericPolicy, + context: StatisticalErrorContext, +) -> Result, StatisticalError> { + if value.is_nan() { + return match policy.missing_data() { + MissingDataPolicy::Propagate => Ok(Some(value)), + MissingDataPolicy::Omit => Ok(None), // * Collection helpers interpret None as “do not include this value” + MissingDataPolicy::Reject => Err(StatisticalError::NonFiniteValue { + context, + value: NonFiniteValue::Nan, + }), + }; + } + + if value.is_infinite() && policy.infinity() == InfinityPolicy::Reject { + return Err(StatisticalError::NonFiniteValue { + context, + value: if value > A::zero() { + NonFiniteValue::PositiveInfinity + } else { + NonFiniteValue::NegativeInfinity + }, + }); + } + + Ok(Some(value)) +} + +/// Applies the shared policy to a pair of values, omitting the pair when one +/// of its values is omitted as missing. +pub(crate) fn apply_numeric_policy_pair( + first: A, + second: A, + policy: NumericPolicy, + context: StatisticalErrorContext, +) -> Result, StatisticalError> { + let first = apply_numeric_policy(first, policy, context)?; + let second = apply_numeric_policy(second, policy, context)?; + match (first, second) { + (Some(first), Some(second)) => Ok(Some((first, second))), + _ => Ok(None), + } +} + +/// Collects values after applying the shared policy while retaining the +/// original iterator index in any error context. +pub(crate) fn collect_with_numeric_policy( + values: I, + policy: NumericPolicy, + context: StatisticalErrorContext, +) -> Result, StatisticalError> +where + A: Float, + I: IntoIterator, +{ + values + .into_iter() + .enumerate() + .try_fold(Vec::new(), |mut collected, (index, value)| { + if let Some(value) = apply_numeric_policy(value, policy, context.with_index(index))? { + collected.push(value); + } + Ok(collected) + }) +} diff --git a/src/summary_statistics/descriptive.rs b/src/summary_statistics/descriptive.rs index 5c87ff2..c1ea896 100644 --- a/src/summary_statistics/descriptive.rs +++ b/src/summary_statistics/descriptive.rs @@ -1,6 +1,7 @@ -use crate::errors::NonFiniteValue; -use crate::errors::SummaryStatisticsError; -use crate::policies::{InfinityPolicy, MissingDataPolicy, NumericPolicy}; +use crate::errors::{ + StatisticalError, StatisticalErrorContext, StatisticalOperation, SummaryStatisticsError, +}; +use crate::policies::{apply_numeric_policy, NumericPolicy}; use num_traits::{Float, FromPrimitive}; use std::cmp::Ordering; @@ -72,47 +73,39 @@ where self.sample_variance().sqrt() } - pub(super) fn from_iter(values: I) -> Result + pub(super) fn from_iter( + values: I, + policy: NumericPolicy, + ) -> Result where I: IntoIterator, { - Self::from_iter_with_policy(values, NumericPolicy::propagate()) + Self::from_iter_with_context( + values, + policy, + StatisticalErrorContext::new(StatisticalOperation::DescriptiveStatistics), + ) + .map_err(legacy_summary_error) } - pub(super) fn from_iter_with_policy( + pub(super) fn from_iter_with_context( values: I, policy: NumericPolicy, - ) -> Result + context: StatisticalErrorContext, + ) -> Result where I: IntoIterator, { let mut accumulator: Option> = None; for (index, value) in values.into_iter().enumerate() { - if value.is_nan() { - match policy.missing_data() { - MissingDataPolicy::Propagate => {} - MissingDataPolicy::Omit => continue, - MissingDataPolicy::Reject => { - return Err(SummaryStatisticsError::NonFiniteValue { - index, - value: NonFiniteValue::Nan, - }); - } - } - } else if value.is_infinite() && policy.infinity() == InfinityPolicy::Reject { - return Err(SummaryStatisticsError::NonFiniteValue { - index, - value: if value > A::zero() { - NonFiniteValue::PositiveInfinity - } else { - NonFiniteValue::NegativeInfinity - }, - }); - } + let Some(value) = apply_numeric_policy(value, policy, context.with_index(index))? + else { + continue; + }; if let Some(ref mut accumulator) = accumulator { - accumulator.update(value)?; + accumulator.update(value, context.with_index(index))?; } else { accumulator = Some(Accumulator::new(value)); } @@ -120,7 +113,7 @@ where accumulator .map(Accumulator::finish) - .ok_or(SummaryStatisticsError::EmptyInput) + .ok_or(StatisticalError::EmptyInput { context }) } fn from_usize(value: usize) -> A { @@ -128,6 +121,24 @@ where } } +fn legacy_summary_error(error: StatisticalError) -> SummaryStatisticsError { + match error { + StatisticalError::EmptyInput { .. } => SummaryStatisticsError::EmptyInput, + StatisticalError::NonFiniteValue { context, value } => { + SummaryStatisticsError::NonFiniteValue { + index: context.index.unwrap_or(0), + value, + } + } + StatisticalError::UndefinedOrder { .. } => SummaryStatisticsError::UndefinedOrder, + StatisticalError::ShapeMismatch { .. } + | StatisticalError::InvalidDegreesOfFreedom { .. } + | StatisticalError::InvalidDomain { .. } => { + unreachable!("descriptive summary cannot produce this error") + } + } +} + struct Accumulator { count: usize, mean: A, @@ -150,7 +161,11 @@ where } } - fn update(&mut self, value: A) -> Result<(), SummaryStatisticsError> { + fn update( + &mut self, + value: A, + context: StatisticalErrorContext, + ) -> Result<(), StatisticalError> { self.count += 1; let count = DescriptiveStatistics::::from_usize(self.count); let delta = value - self.mean; @@ -160,12 +175,12 @@ where match value.partial_cmp(&self.min) { Some(Ordering::Less) => self.min = value, Some(_) => {} - None => return Err(SummaryStatisticsError::UndefinedOrder), + None => return Err(StatisticalError::UndefinedOrder { context }), } match value.partial_cmp(&self.max) { Some(Ordering::Greater) => self.max = value, Some(_) => {} - None => return Err(SummaryStatisticsError::UndefinedOrder), + None => return Err(StatisticalError::UndefinedOrder { context }), } Ok(()) diff --git a/src/summary_statistics/means.rs b/src/summary_statistics/means.rs index af47e43..03ebead 100644 --- a/src/summary_statistics/means.rs +++ b/src/summary_statistics/means.rs @@ -1,6 +1,10 @@ use super::DescriptiveStatistics; use super::SummaryStatisticsExt; -use crate::errors::{EmptyInput, MultiInputError, ShapeMismatch, SummaryStatisticsError}; +use crate::errors::{ + EmptyInput, MultiInputError, ShapeMismatch, StatisticalError, StatisticalErrorContext, + StatisticalOperation, SummaryStatisticsError, +}; +use crate::policies::{apply_numeric_policy_pair, collect_with_numeric_policy, NumericPolicy}; use ndarray::{Array, ArrayBase, ArrayRef, Axis, Data, Dimension, Ix1, RemoveAxis}; use num_integer::IterBinomial; use num_traits::{Float, FromPrimitive, Zero}; @@ -14,7 +18,21 @@ where where A: Float + FromPrimitive, { - DescriptiveStatistics::from_iter(self.iter().copied()) + DescriptiveStatistics::from_iter(self.iter().copied(), NumericPolicy::propagate()) + } + + fn descriptive_statistics_with_policy( + &self, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + { + DescriptiveStatistics::from_iter_with_context( + self.iter().copied(), + policy, + StatisticalErrorContext::new(StatisticalOperation::DescriptiveStatistics), + ) } fn descriptive_statistics_axis( @@ -33,47 +51,363 @@ where let summaries = self .lanes(axis) .into_iter() - .map(|lane| DescriptiveStatistics::from_iter(lane.iter().copied())) + .map(|lane| { + DescriptiveStatistics::from_iter(lane.iter().copied(), NumericPolicy::propagate()) + }) .collect::, _>>()?; Ok(Array::from_shape_vec(shape, summaries) .expect("descriptive-statistics lanes must match the output shape")) } - fn descriptive_statistics_with_policy( - &self, - policy: crate::policies::NumericPolicy, - ) -> Result, SummaryStatisticsError> - where - A: Float + FromPrimitive, - { - DescriptiveStatistics::from_iter_with_policy(self.iter().copied(), policy) - } - fn descriptive_statistics_axis_with_policy( &self, - axis: ndarray::Axis, - policy: crate::policies::NumericPolicy, - ) -> Result, D::Smaller>, SummaryStatisticsError> + axis: Axis, + policy: NumericPolicy, + ) -> Result, D::Smaller>, StatisticalError> where A: Float + FromPrimitive, D: RemoveAxis, { if self.is_empty() { - return Err(SummaryStatisticsError::EmptyInput); + return Err(StatisticalError::EmptyInput { + context: StatisticalErrorContext::axis( + StatisticalOperation::DescriptiveStatistics, + axis.index(), + ), + }); } let shape = self.raw_dim().remove_axis(axis); let summaries = self .lanes(axis) .into_iter() - .map(|lane| DescriptiveStatistics::from_iter_with_policy(lane.iter().copied(), policy)) + .enumerate() + .map(|(lane_index, lane)| { + DescriptiveStatistics::from_iter_with_context( + lane.iter().copied(), + policy, + StatisticalErrorContext::axis_lane( + StatisticalOperation::DescriptiveStatistics, + axis.index(), + lane_index, + ), + ) + }) .collect::, _>>()?; Ok(Array::from_shape_vec(shape, summaries) .expect("descriptive-statistics lanes must match the output shape")) } + fn mean_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::Mean); + let values = collect_with_numeric_policy(self.iter().copied(), policy, context)?; + mean_values(&values, context) + } + + fn weighted_mean_with_policy( + &self, + weights: &Self, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::WeightedMean); + let pairs = collect_weight_pairs(self, weights, policy, context)?; + weighted_mean_values(&pairs, context) + } + + fn weighted_sum_with_policy( + &self, + weights: &Self, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::WeightedSum); + let pairs = collect_weight_pairs(self, weights, policy, context)?; + weighted_sum_values(&pairs, context) + } + + fn weighted_mean_axis_with_policy( + &self, + axis: Axis, + weights: &ArrayRef, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + D: RemoveAxis, + { + let context = + StatisticalErrorContext::axis(StatisticalOperation::WeightedMean, axis.index()); + let pairs = weighted_axis_values(self, axis, weights, policy, context)?; + let shape = self.raw_dim().remove_axis(axis); + let values = pairs + .into_iter() + .enumerate() + .map(|(lane, pairs)| { + weighted_mean_values( + &pairs, + StatisticalErrorContext::axis_lane( + StatisticalOperation::WeightedMean, + axis.index(), + lane, + ), + ) + }) + .collect::, _>>()?; + Ok(Array::from_shape_vec(shape, values) + .expect("weighted mean lanes must match the output shape")) + } + + fn weighted_sum_axis_with_policy( + &self, + axis: Axis, + weights: &ArrayRef, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + D: RemoveAxis, + { + let context = + StatisticalErrorContext::axis(StatisticalOperation::WeightedSum, axis.index()); + let pairs = weighted_axis_values(self, axis, weights, policy, context)?; + let shape = self.raw_dim().remove_axis(axis); + let values = pairs + .into_iter() + .enumerate() + .map(|(lane, pairs)| { + weighted_sum_values( + &pairs, + StatisticalErrorContext::axis_lane( + StatisticalOperation::WeightedSum, + axis.index(), + lane, + ), + ) + }) + .collect::, _>>()?; + Ok(Array::from_shape_vec(shape, values) + .expect("weighted sum lanes must match the output shape")) + } + + fn harmonic_mean_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::HarmonicMean); + let values = collect_with_numeric_policy(self.iter().copied(), policy, context)?; + if values.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + let sum = values + .iter() + .copied() + .fold(A::zero(), |sum, value| sum + value.recip()); + Ok((A::from_usize(values.len()).expect("length must convert to A") / sum).recip()) + } + + fn geometric_mean_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::GeometricMean); + let values = collect_with_numeric_policy(self.iter().copied(), policy, context)?; + if values.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + let sum = values + .iter() + .copied() + .fold(A::zero(), |sum, value| sum + value.ln()); + Ok((sum / A::from_usize(values.len()).expect("length must convert to A")).exp()) + } + + fn raw_moment_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive, + { + Ok(self.raw_moments_with_policy(order, policy)?[usize::from(order)]) + } + + fn raw_moments_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::RawMoment); + let values = collect_with_numeric_policy(self.iter().copied(), policy, context)?; + moments_values(&values, order, context) + } + + fn standardized_moment_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive, + { + Ok(self.standardized_moments_with_policy(order, policy)?[usize::from(order)]) + } + + fn standardized_moments_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + { + let central_moments = self.central_moments_with_policy(order, policy)?; + if order < 2 { + return Ok(central_moments); + } + let standard_deviation = central_moments[2].sqrt(); + Ok(central_moments + .into_iter() + .enumerate() + .map(|(order, moment)| { + if order < 2 { + moment + } else { + moment / standard_deviation.powi(order as i32) + } + }) + .collect()) + } + + fn central_moment_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive, + { + Ok(self.central_moments_with_policy(order, policy)?[usize::from(order)]) + } + + fn central_moments_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::CentralMoment); + let values = collect_with_numeric_policy(self.iter().copied(), policy, context)?; + central_moments_values(&values, order, context) + } + + fn weighted_var_with_policy( + &self, + weights: &Self, + ddof: A, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::WeightedVariance); + let pairs = collect_weight_pairs(self, weights, policy, context)?; + weighted_variance_values(&pairs, ddof, context) + } + + fn weighted_std_with_policy( + &self, + weights: &Self, + ddof: A, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive, + { + Ok(self.weighted_var_with_policy(weights, ddof, policy)?.sqrt()) + } + + fn weighted_var_axis_with_policy( + &self, + axis: Axis, + weights: &ArrayRef, + ddof: A, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + D: RemoveAxis, + { + let context = + StatisticalErrorContext::axis(StatisticalOperation::WeightedVariance, axis.index()); + let pairs = weighted_axis_values(self, axis, weights, policy, context)?; + let shape = self.raw_dim().remove_axis(axis); + let values = pairs + .into_iter() + .enumerate() + .map(|(lane, pairs)| { + weighted_variance_values( + &pairs, + ddof, + StatisticalErrorContext::axis_lane( + StatisticalOperation::WeightedVariance, + axis.index(), + lane, + ), + ) + }) + .collect::, _>>()?; + Ok(Array::from_shape_vec(shape, values) + .expect("weighted variance lanes must match the output shape")) + } + + fn weighted_std_axis_with_policy( + &self, + axis: Axis, + weights: &ArrayRef, + ddof: A, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + D: RemoveAxis, + { + Ok(self + .weighted_var_axis_with_policy(axis, weights, ddof, policy)? + .mapv_into(|value| value.sqrt())) + } + + fn kurtosis_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive, + { + let moments = self.central_moments_with_policy(4, policy)?; + Ok(moments[4] / moments[2].powi(2)) + } + + fn skewness_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive, + { + let moments = self.central_moments_with_policy(3, policy)?; + Ok(moments[3] / moments[2].sqrt().powi(3)) + } + fn mean(&self) -> Result where A: Clone + FromPrimitive + Add + Div + Zero, @@ -391,6 +725,226 @@ where private_impl! {} } +fn mean_values(values: &[A], context: StatisticalErrorContext) -> Result +where + A: Float + FromPrimitive, +{ + if values.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + let sum = values + .iter() + .copied() + .fold(A::zero(), |sum, value| sum + value); + Ok(sum / A::from_usize(values.len()).expect("length must convert to A")) +} + +fn collect_weight_pairs( + values: &ArrayRef, + weights: &ArrayRef, + policy: NumericPolicy, + context: StatisticalErrorContext, +) -> Result, StatisticalError> +where + A: Float, + D: Dimension, +{ + if values.shape() != weights.shape() { + return Err(StatisticalError::ShapeMismatch { + context, + first_shape: values.shape().to_vec(), + second_shape: weights.shape().to_vec(), + }); + } + + values + .iter() + .copied() + .zip(weights.iter().copied()) + .enumerate() + .try_fold(Vec::new(), |mut pairs, (index, (value, weight))| { + if let Some(pair) = + apply_numeric_policy_pair(value, weight, policy, context.with_index(index))? + { + pairs.push(pair); + } + Ok(pairs) + }) +} + +fn weighted_axis_values( + values: &ArrayRef, + axis: Axis, + weights: &ArrayRef, + policy: NumericPolicy, + context: StatisticalErrorContext, +) -> Result>, StatisticalError> +where + A: Float, + D: Dimension, +{ + if values.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + if values.shape()[axis.index()] != weights.len() { + return Err(StatisticalError::ShapeMismatch { + context, + first_shape: values.shape().to_vec(), + second_shape: weights.shape().to_vec(), + }); + } + + values + .lanes(axis) + .into_iter() + .enumerate() + .map(|(lane, values)| { + values + .iter() + .copied() + .zip(weights.iter().copied()) + .enumerate() + .try_fold(Vec::new(), |mut pairs, (index, (value, weight))| { + if let Some(pair) = apply_numeric_policy_pair( + value, + weight, + policy, + StatisticalErrorContext::axis_lane(context.operation, axis.index(), lane) + .with_index(index), + )? { + pairs.push(pair); + } + Ok(pairs) + }) + }) + .collect() +} + +fn weighted_sum_values( + pairs: &[(A, A)], + context: StatisticalErrorContext, +) -> Result +where + A: Float, +{ + if pairs.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(pairs + .iter() + .copied() + .fold(A::zero(), |sum, (value, weight)| sum + value * weight)) +} + +fn weighted_mean_values( + pairs: &[(A, A)], + context: StatisticalErrorContext, +) -> Result +where + A: Float, +{ + let weighted_sum = weighted_sum_values(pairs, context)?; + let weight_sum = pairs + .iter() + .copied() + .fold(A::zero(), |sum, (_, weight)| sum + weight); + Ok(weighted_sum / weight_sum) +} + +fn moments_values( + values: &[A], + order: u16, + context: StatisticalErrorContext, +) -> Result, StatisticalError> +where + A: Float + FromPrimitive, +{ + if values.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + + let n = A::from_usize(values.len()).expect("length must convert to A"); + let order = i32::from(order); + let mut moments = vec![A::one()]; + if order >= 1 { + moments.push( + values + .iter() + .copied() + .fold(A::zero(), |sum, value| sum + value) + / n, + ); + } + for power in 2..=order { + moments.push( + values + .iter() + .copied() + .fold(A::zero(), |sum, value| sum + value.powi(power)) + / n, + ); + } + Ok(moments) +} + +fn central_moments_values( + values: &[A], + order: u16, + context: StatisticalErrorContext, +) -> Result, StatisticalError> +where + A: Float + FromPrimitive, +{ + if values.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + match order { + 0 => Ok(vec![A::one()]), + 1 => Ok(vec![A::one(), A::zero()]), + n => { + let mean = mean_values(values, context)?; + let shifted: Vec = values.iter().copied().map(|value| value - mean).collect(); + let shifted_moments = moments_values(&shifted, n, context)?; + let correction_term = -shifted_moments[1]; + let mut central_moments = vec![A::one(), A::zero()]; + for k in 2..=n { + let coefficients = central_moment_coefficients(&shifted_moments[..=(k as usize)]); + central_moments.push(horner_method(coefficients, correction_term)); + } + Ok(central_moments) + } + } +} + +fn weighted_variance_values( + pairs: &[(A, A)], + ddof: A, + context: StatisticalErrorContext, +) -> Result +where + A: Float + FromPrimitive, +{ + if pairs.is_empty() { + return Err(StatisticalError::EmptyInput { context }); + } + let zero = A::zero(); + let one = A::one(); + if ddof < zero || ddof > one { + return Err(StatisticalError::InvalidDegreesOfFreedom { context }); + } + + let mut weight_sum = zero; + let mut mean = zero; + let mut sum = zero; + for &(value, weight) in pairs { + weight_sum = weight_sum + weight; + let value_minus_mean = value - mean; + mean = mean + (weight / weight_sum) * value_minus_mean; + sum = sum + weight * value_minus_mean * (value - mean); + } + Ok(sum / (weight_sum - ddof)) +} + /// Private function for `weighted_var` without conditions and asserts. fn inner_weighted_var( arr: &ArrayRef, diff --git a/src/summary_statistics/mod.rs b/src/summary_statistics/mod.rs index 7c6b05b..ee15b1c 100644 --- a/src/summary_statistics/mod.rs +++ b/src/summary_statistics/mod.rs @@ -1,5 +1,5 @@ //! Summary statistics (e.g. mean, variance, etc.). -use crate::errors::{EmptyInput, MultiInputError, SummaryStatisticsError}; +use crate::errors::{EmptyInput, MultiInputError, StatisticalError, SummaryStatisticsError}; use crate::policies::NumericPolicy; use ndarray::{Array, ArrayRef, Axis, Dimension, Ix1, RemoveAxis}; use num_traits::{Float, FromPrimitive, Zero}; @@ -15,43 +15,17 @@ pub trait SummaryStatisticsExt where D: Dimension, { - /// Returns a fused descriptive-statistics summary of all elements in the array. + /// Returns a fused descriptive-statistics summary of all elements in the + /// array using the historical propagate behavior. /// - /// If the array is empty, `SummaryStatisticsError::EmptyInput` is returned. - /// If a required minimum or maximum comparison has undefined ordering, - /// `SummaryStatisticsError::UndefinedOrder` is returned. - fn descriptive_statistics(&self) -> Result, SummaryStatisticsError> - where - A: Float + FromPrimitive; - - /// Returns a descriptive-statistics summary for every lane along `axis`. - /// - /// The returned array has the input shape with `axis` removed. The method - /// panics if `axis` is out of bounds and returns the first summary error - /// encountered while processing the lanes. - fn descriptive_statistics_axis( - &self, - axis: Axis, - ) -> Result, D::Smaller>, SummaryStatisticsError> - where - A: Float + FromPrimitive, - D: RemoveAxis; - - /// Returns a fused descriptive-statistics summary using an explicit - /// missing-data and infinity policy. - /// - /// [`NumericPolicy::propagate`] preserves the behavior of - /// [`descriptive_statistics`](Self::descriptive_statistics). Use - /// [`NumericPolicy::omit_missing`] to omit NaN values. Infinity is not - /// treated as missing: it is propagated by that policy and can instead be - /// rejected with [`NumericPolicy::reject_non_finite`]. If omission removes - /// every value, `SummaryStatisticsError::EmptyInput` is returned. + /// Use [`Self::descriptive_statistics_with_policy`] with + /// [`NumericPolicy::omit_missing`] for explicit omission semantics. /// /// # Example /// /// ``` /// use ndarray::array; - /// use ndarray_stats::{NumericPolicy, SummaryStatisticsExt}; + /// use ndarray_stats::{policies::NumericPolicy, SummaryStatisticsExt}; /// /// let summary = array![1.0, f64::NAN, 3.0] /// .descriptive_statistics_with_policy(NumericPolicy::omit_missing()) @@ -59,25 +33,56 @@ where /// assert_eq!(summary.count(), 2); /// assert_eq!(summary.mean(), 2.0); /// ``` + fn descriptive_statistics(&self) -> Result, SummaryStatisticsError> + where + A: Float + FromPrimitive; + + /// Returns a fused descriptive-statistics summary under an explicit + /// [`NumericPolicy`]. Policy-aware methods use the shared + /// [`StatisticalError`] model. fn descriptive_statistics_with_policy( &self, policy: NumericPolicy, - ) -> Result, SummaryStatisticsError> + ) -> Result, StatisticalError> where A: Float + FromPrimitive; /// Returns a descriptive-statistics summary for every lane along `axis` - /// using an explicit missing-data and infinity policy. + /// using the historical propagate behavior. + /// + /// The returned array has the input shape with `axis` removed. The method + /// panics if `axis` is out of bounds and returns the first summary error + /// Use [`Self::descriptive_statistics_axis_with_policy`] for explicit + /// omission or rejection semantics. + /// + /// # Example + /// + /// ``` + /// use ndarray::{array, Axis}; + /// use ndarray_stats::{policies::NumericPolicy, SummaryStatisticsExt}; /// - /// The returned array has the input shape with `axis` removed. Omission is - /// applied independently to each lane. If any lane is empty after - /// omission, `SummaryStatisticsError::EmptyInput` is returned because the - /// current summary result type cannot represent an absent lane. + /// let summaries = array![[1.0, f64::NAN, 3.0], [4.0, 5.0, 6.0]] + /// .descriptive_statistics_axis_with_policy(Axis(1), NumericPolicy::omit_missing()) + /// .unwrap(); + /// assert_eq!(summaries[0].count(), 2); + /// assert_eq!(summaries[0].mean(), 2.0); + /// ``` + fn descriptive_statistics_axis( + &self, + axis: Axis, + ) -> Result, D::Smaller>, SummaryStatisticsError> + where + A: Float + FromPrimitive, + D: RemoveAxis; + + /// Returns summaries for every lane along `axis` under an explicit + /// [`NumericPolicy`]. An omitted lane still returns the legacy empty-input + /// error; a per-lane diagnostic result is deferred as an API decision. fn descriptive_statistics_axis_with_policy( &self, axis: Axis, policy: NumericPolicy, - ) -> Result, D::Smaller>, SummaryStatisticsError> + ) -> Result, D::Smaller>, StatisticalError> where A: Float + FromPrimitive, D: RemoveAxis; @@ -119,7 +124,7 @@ where /// * `MultiInputError::EmptyInput` if `self` is empty /// * `MultiInputError::ShapeMismatch` if `self` and `weights` don't have the same shape /// - /// [`arithmetic weighted mean`] https://en.wikipedia.org/wiki/Weighted_arithmetic_mean + /// [`arithmetic weighted mean`] fn weighted_mean(&self, weights: &Self) -> Result where A: Copy + Div + Mul + Zero; @@ -160,7 +165,7 @@ where /// * `MultiInputError::EmptyInput` if `self` is empty /// * `MultiInputError::ShapeMismatch` if `self` length along axis is not equal to `weights` length /// - /// [`arithmetic weighted mean`] https://en.wikipedia.org/wiki/Weighted_arithmetic_mean + /// [`arithmetic weighted mean`] fn weighted_mean_axis( &self, axis: Axis, @@ -466,6 +471,154 @@ where where A: Float + FromPrimitive; + /// Policy-aware arithmetic mean. + fn mean_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive; + + /// Policy-aware weighted mean and weighted sum. + fn weighted_mean_with_policy( + &self, + weights: &Self, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive; + + fn weighted_sum_with_policy( + &self, + weights: &Self, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive; + + /// Policy-aware weighted reductions along an axis. + fn weighted_mean_axis_with_policy( + &self, + axis: Axis, + weights: &ArrayRef, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + D: RemoveAxis; + + fn weighted_sum_axis_with_policy( + &self, + axis: Axis, + weights: &ArrayRef, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + D: RemoveAxis; + + /// Policy-aware scalar summaries and moments. + fn harmonic_mean_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive; + + fn geometric_mean_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive; + + fn raw_moment_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive; + + fn raw_moments_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive; + + fn standardized_moment_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive; + + fn standardized_moments_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive; + + fn central_moment_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive; + + fn central_moments_with_policy( + &self, + order: u16, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive; + + fn weighted_var_with_policy( + &self, + weights: &Self, + ddof: A, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive; + + fn weighted_std_with_policy( + &self, + weights: &Self, + ddof: A, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive; + + fn weighted_var_axis_with_policy( + &self, + axis: Axis, + weights: &ArrayRef, + ddof: A, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + D: RemoveAxis; + + fn weighted_std_axis_with_policy( + &self, + axis: Axis, + weights: &ArrayRef, + ddof: A, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + D: RemoveAxis; + + fn kurtosis_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive; + + fn skewness_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive; + private_decl! {} } diff --git a/tests/numeric_policy.rs b/tests/numeric_policy.rs new file mode 100644 index 0000000..fba410d --- /dev/null +++ b/tests/numeric_policy.rs @@ -0,0 +1,103 @@ +use approx::assert_abs_diff_eq; +use ndarray::{arr2, array}; +use ndarray_stats::{ + errors::{StatisticalError, StatisticalOperation}, + policies::NumericPolicy, + CorrelationExt, DeviationExt, EntropyExt, SummaryStatisticsExt, +}; + +#[test] +fn omission_is_shared_by_means_moments_and_weighted_reductions() { + let data = array![1.0, f64::NAN, 3.0]; + let weights = array![1.0, f64::NAN, 1.0]; + let policy = NumericPolicy::omit_missing(); + + assert_eq!(data.mean_with_policy(policy).unwrap(), 2.0); + assert_eq!(data.raw_moment_with_policy(2, policy).unwrap(), 5.0); + assert_eq!(data.central_moment_with_policy(2, policy).unwrap(), 1.0); + assert_eq!( + data.weighted_sum_with_policy(&weights, policy).unwrap(), + 4.0 + ); + assert_eq!( + data.weighted_mean_with_policy(&weights, policy).unwrap(), + 2.0 + ); +} + +#[test] +fn omission_is_shared_by_entropy_and_deviations() { + let data = array![0.5, f64::NAN, 0.5]; + let other = array![0.5, 0.25, 0.25]; + let policy = NumericPolicy::omit_missing(); + + assert_abs_diff_eq!(data.entropy_with_policy(policy).unwrap(), 2f64.ln()); + assert_eq!( + array![1.0, f64::NAN, 3.0] + .count_eq_with_policy(&array![1.0, 2.0, 2.0], policy) + .unwrap(), + 1 + ); + assert_eq!( + array![1.0, f64::NAN, 3.0] + .sq_l2_dist_with_policy(&array![1.0, 2.0, 2.0], policy) + .unwrap(), + 1.0 + ); + assert!(data.kl_divergence_with_policy(&other, policy).is_ok()); +} + +#[test] +fn covariance_uses_listwise_policy_and_preserves_correlation() { + let data = arr2(&[[1.0, f64::NAN, 3.0], [1.0, 2.0, 3.0]]); + let policy = NumericPolicy::omit_missing(); + + let covariance = data.cov_with_policy(0.0, policy).unwrap(); + assert_abs_diff_eq!(covariance[[0, 0]], 1.0, epsilon = 1e-12); + assert_abs_diff_eq!(covariance[[0, 1]], 1.0, epsilon = 1e-12); + assert_abs_diff_eq!( + data.pearson_correlation_with_policy(policy).unwrap()[[0, 1]], + 1.0, + epsilon = 1e-12 + ); +} + +#[test] +fn rejection_reports_operation_lane_and_original_index() { + let data = arr2(&[[1.0, f64::NAN], [2.0, 3.0]]); + let error = data + .cov_with_policy(0.0, NumericPolicy::reject_non_finite()) + .unwrap_err(); + + assert_eq!( + error, + StatisticalError::NonFiniteValue { + context: crate_context(StatisticalOperation::Covariance, 1, 0, 1), + value: ndarray_stats::errors::NonFiniteValue::Nan, + } + ); +} + +#[test] +fn infinity_remains_distinct_from_missing_values() { + let error = array![1.0, f64::INFINITY] + .mean_with_policy(NumericPolicy::reject_non_finite()) + .unwrap_err(); + assert!(matches!( + error, + StatisticalError::NonFiniteValue { + context, + value: ndarray_stats::errors::NonFiniteValue::PositiveInfinity, + } if context.operation == StatisticalOperation::Mean && context.index == Some(1) + )); +} + +fn crate_context( + operation: StatisticalOperation, + axis: usize, + lane: usize, + index: usize, +) -> ndarray_stats::errors::StatisticalErrorContext { + ndarray_stats::errors::StatisticalErrorContext::axis_lane(operation, axis, lane) + .with_index(index) +} diff --git a/tests/policies.rs b/tests/policies.rs index e646a67..0f696ec 100644 --- a/tests/policies.rs +++ b/tests/policies.rs @@ -1,8 +1,9 @@ use approx::assert_abs_diff_eq; use ndarray::{array, Axis}; use ndarray_stats::{ - errors::{NonFiniteValue, SummaryStatisticsError}, - InfinityPolicy, MissingDataPolicy, NumericPolicy, SummaryStatisticsExt, + errors::{NonFiniteValue, StatisticalError, StatisticalErrorContext, StatisticalOperation}, + policies::{InfinityPolicy, MissingDataPolicy, NumericPolicy}, + SummaryStatisticsExt, }; #[test] @@ -41,7 +42,9 @@ fn omission_of_every_value_is_an_empty_input() { assert_eq!( data.descriptive_statistics_with_policy(NumericPolicy::omit_missing()), - Err(SummaryStatisticsError::EmptyInput) + Err(StatisticalError::EmptyInput { + context: StatisticalErrorContext::new(StatisticalOperation::DescriptiveStatistics), + }) ); } @@ -51,7 +54,13 @@ fn omission_of_every_value_in_one_axis_lane_is_an_empty_input() { assert_eq!( data.descriptive_statistics_axis_with_policy(Axis(1), NumericPolicy::omit_missing()), - Err(SummaryStatisticsError::EmptyInput) + Err(StatisticalError::EmptyInput { + context: StatisticalErrorContext::axis_lane( + StatisticalOperation::DescriptiveStatistics, + 1, + 1, + ), + }) ); } @@ -77,8 +86,9 @@ fn reject_policy_reports_nan_and_infinity() { let nan_result = array![1.0, f64::NAN].descriptive_statistics_with_policy(reject_nan); assert_eq!( nan_result, - Err(SummaryStatisticsError::NonFiniteValue { - index: 1, + Err(StatisticalError::NonFiniteValue { + context: StatisticalErrorContext::new(StatisticalOperation::DescriptiveStatistics) + .with_index(1), value: NonFiniteValue::Nan, }) ); @@ -87,8 +97,9 @@ fn reject_policy_reports_nan_and_infinity() { .descriptive_statistics_with_policy(NumericPolicy::reject_non_finite()); assert_eq!( infinity_result, - Err(SummaryStatisticsError::NonFiniteValue { - index: 1, + Err(StatisticalError::NonFiniteValue { + context: StatisticalErrorContext::new(StatisticalOperation::DescriptiveStatistics) + .with_index(1), value: NonFiniteValue::NegativeInfinity, }) ); @@ -100,6 +111,9 @@ fn default_policy_preserves_existing_nan_behavior() { assert_eq!( data.descriptive_statistics_with_policy(NumericPolicy::default()), - data.descriptive_statistics() + Err(StatisticalError::UndefinedOrder { + context: StatisticalErrorContext::new(StatisticalOperation::DescriptiveStatistics) + .with_index(1), + }) ); } From 7477e0c3257e908535eb69f990ad65b3119a1054 Mon Sep 17 00:00:00 2001 From: N3uralN3twork Date: Sun, 2 Aug 2026 19:38:57 -0400 Subject: [PATCH 2/2] Fixing the PR gaps --- benches/summary_statistics.rs | 31 +++ docs/SUMMARY_STATISTICS.md | 8 +- src/correlation.rs | 28 ++- src/deviation.rs | 175 ++++++++++------ src/entropy.rs | 77 +++---- src/policies.rs | 48 ++++- src/summary_statistics/means.rs | 357 ++++++++++++++++++++------------ tests/numeric_policy.rs | 205 +++++++++++++++++- 8 files changed, 674 insertions(+), 255 deletions(-) diff --git a/benches/summary_statistics.rs b/benches/summary_statistics.rs index 1f9013a..c8786d3 100644 --- a/benches/summary_statistics.rs +++ b/benches/summary_statistics.rs @@ -78,6 +78,13 @@ fn descriptive_statistics(c: &mut Criterion) { let weights = Array1::ones(len); group.bench_with_input(format!("fused/{len}"), &data, |b, data| { + b.iter(|| { + let summary = black_box(data.descriptive_statistics().unwrap()); + black_box(score_f64(&summary)); + }) + }); + + group.bench_with_input(format!("policy/{len}"), &data, |b, data| { b.iter(|| { let summary = black_box( data.descriptive_statistics_with_policy(NumericPolicy::default()) @@ -117,6 +124,14 @@ fn descriptive_statistics_axis(c: &mut Criterion) { let axis_one_weights = Array1::ones(8); group.bench_with_input(format!("fused_axis0/{len}"), &data, |b, data| { + b.iter(|| { + let summaries = black_box(data.descriptive_statistics_axis(Axis(0)).unwrap()); + let score = summaries.iter().map(score_f64).sum::(); + black_box(score); + }) + }); + + group.bench_with_input(format!("policy_axis0/{len}"), &data, |b, data| { b.iter(|| { let summaries = black_box( data.descriptive_statistics_axis_with_policy(Axis(0), NumericPolicy::default()) @@ -134,6 +149,14 @@ fn descriptive_statistics_axis(c: &mut Criterion) { }); group.bench_with_input(format!("fused_axis1/{len}"), &data, |b, data| { + b.iter(|| { + let summaries = black_box(data.descriptive_statistics_axis(Axis(1)).unwrap()); + let score = summaries.iter().map(score_f64).sum::(); + black_box(score); + }) + }); + + group.bench_with_input(format!("policy_axis1/{len}"), &data, |b, data| { b.iter(|| { let summaries = black_box( data.descriptive_statistics_axis_with_policy(Axis(1), NumericPolicy::default()) @@ -161,6 +184,14 @@ fn descriptive_statistics_f32(c: &mut Criterion) { for &len in &lens { group.bench_function(format!("fused/{len}"), |b| { + let data: Array1 = Array::random(len, Uniform::new(0.0, 1.0).unwrap()); + b.iter(|| { + let summary = black_box(data.descriptive_statistics().unwrap()); + black_box(score_f32(&summary)); + }) + }); + + group.bench_function(format!("policy/{len}"), |b| { let data: Array1 = Array::random(len, Uniform::new(0.0, 1.0).unwrap()); b.iter(|| { let summary = black_box( diff --git a/docs/SUMMARY_STATISTICS.md b/docs/SUMMARY_STATISTICS.md index 8b69789..22d75e1 100644 --- a/docs/SUMMARY_STATISTICS.md +++ b/docs/SUMMARY_STATISTICS.md @@ -175,11 +175,9 @@ With `NumericPolicy::reject_non_finite()`, the same calls return a `NumericPolicy::propagate()`, no filtering occurs and IEEE-754 behavior is preserved. -The current `harmonic_mean_with_policy` implementation returns the reciprocal -of the intended harmonic mean (`0.666...` for `missing`), while the legacy -`harmonic_mean` returns `1.5`. The table documents the method's stated -statistical meaning; this implementation discrepancy is not changed by this -documentation-only update. +`harmonic_mean_with_policy` applies the selected policy before calculating +`n / sum(1 / x)`, so the `missing` example above returns `1.5`, matching the +legacy `harmonic_mean` result for the retained values. ## Errors and edge cases diff --git a/src/correlation.rs b/src/correlation.rs index f89d22b..f842715 100644 --- a/src/correlation.rs +++ b/src/correlation.rs @@ -365,38 +365,36 @@ where return Ok(Array2::zeros((0, n_observations))); } - let mut filtered = vec![Vec::new(); n_variables]; + if policy == NumericPolicy::propagate() { + return Ok(data.to_owned()); + } + + let mut retained_observations = Vec::with_capacity(n_observations); for observation in 0..n_observations { - let mut values = Vec::with_capacity(n_variables); let mut include = true; for variable in 0..n_variables { let context = StatisticalErrorContext::axis_lane(operation, 1, variable).with_index(observation); match apply_numeric_policy(data[[variable, observation]], policy, context)? { - Some(value) => values.push(value), + Some(_) => {} None => { include = false; - values.push(A::zero()); } } } if include { - for variable in 0..n_variables { - filtered[variable].push(values[variable]); - } + retained_observations.push(observation); } } - let retained = filtered[0].len(); + let retained = retained_observations.len(); if retained == 0 { return Err(StatisticalError::EmptyInput { context }); } - let values = filtered - .into_iter() - .flat_map(|row| row.into_iter()) - .collect(); - Ok(Array2::from_shape_vec((n_variables, retained), values) - .expect("filtered covariance rows must have equal lengths")) + Ok(Array2::from_shape_fn( + (n_variables, retained), + |(variable, index)| data[[variable, retained_observations[index]]], + )) } fn covariance_with_policy( @@ -412,7 +410,7 @@ where let filtered = policy_matrix(data, policy, operation)?; let (n_variables, n_observations) = filtered.dim(); let n_observations_as_a = A::from_usize(n_observations).expect("length must convert to A"); - if ddof < A::zero() || ddof >= n_observations_as_a { + if ddof.is_nan() || ddof < A::zero() || ddof >= n_observations_as_a { return Err(StatisticalError::InvalidDegreesOfFreedom { context }); } if n_variables == 0 { diff --git a/src/deviation.rs b/src/deviation.rs index 5c30a2a..ea1f2a8 100644 --- a/src/deviation.rs +++ b/src/deviation.rs @@ -445,13 +445,21 @@ where where A: Float, { - let pairs = collect_policy_pairs( + let (count, included) = fold_policy_pairs( self, other, policy, StatisticalErrorContext::new(StatisticalOperation::Deviation), + 0usize, + |count, first, second| count + usize::from(first == second), )?; - Ok(pairs.iter().filter(|(a, b)| a == b).count()) + if included == 0 { + Err(StatisticalError::EmptyInput { + context: StatisticalErrorContext::new(StatisticalOperation::Deviation), + }) + } else { + Ok(count) + } } fn count_neq_with_policy( @@ -462,13 +470,21 @@ where where A: Float, { - let pairs = collect_policy_pairs( + let (count, included) = fold_policy_pairs( self, other, policy, StatisticalErrorContext::new(StatisticalOperation::Deviation), + 0usize, + |count, first, second| count + usize::from(first != second), )?; - Ok(pairs.iter().filter(|(a, b)| a != b).count()) + if included == 0 { + Err(StatisticalError::EmptyInput { + context: StatisticalErrorContext::new(StatisticalOperation::Deviation), + }) + } else { + Ok(count) + } } fn sq_l2_dist_with_policy( @@ -480,14 +496,21 @@ where A: Float, { let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); - let pairs = collect_policy_pairs(self, other, policy, context)?; - if pairs.is_empty() { + let (sum, count) = fold_policy_pairs( + self, + other, + policy, + context, + A::zero(), + |sum, first, second| { + let difference = first - second; + sum + difference * difference + }, + )?; + if count == 0 { return Err(StatisticalError::EmptyInput { context }); } - Ok(pairs.iter().copied().fold(A::zero(), |sum, (a, b)| { - let difference = a - b; - sum + difference * difference - })) + Ok(sum) } fn l2_dist_with_policy( @@ -498,11 +521,22 @@ where where A: Float, { - Ok(self - .sq_l2_dist_with_policy(other, policy)? - .to_f64() - .expect("failed cast from type A to f64") - .sqrt()) + let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); + let (sum, count) = fold_policy_pairs( + self, + other, + policy, + context, + A::zero(), + |sum, first, second| { + let difference = first - second; + sum + difference * difference + }, + )?; + if count == 0 { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(sum.to_f64().expect("failed cast from type A to f64").sqrt()) } fn l1_dist_with_policy( @@ -514,14 +548,18 @@ where A: Float, { let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); - let pairs = collect_policy_pairs(self, other, policy, context)?; - if pairs.is_empty() { + let (sum, count) = fold_policy_pairs( + self, + other, + policy, + context, + A::zero(), + |sum, first, second| sum + (first - second).abs(), + )?; + if count == 0 { return Err(StatisticalError::EmptyInput { context }); } - Ok(pairs - .iter() - .copied() - .fold(A::zero(), |sum, (a, b)| sum + (a - b).abs())) + Ok(sum) } fn linf_dist_with_policy( @@ -533,18 +571,25 @@ where A: Float, { let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); - let pairs = collect_policy_pairs(self, other, policy, context)?; - if pairs.is_empty() { + let (maximum, count) = fold_policy_pairs( + self, + other, + policy, + context, + A::zero(), + |maximum, first, second| { + let value = (first - second).abs(); + if value > maximum { + value + } else { + maximum + } + }, + )?; + if count == 0 { return Err(StatisticalError::EmptyInput { context }); } - Ok(pairs - .iter() - .copied() - .map(|(a, b)| (a - b).abs()) - .fold( - A::zero(), - |max, value| if value > max { value } else { max }, - )) + Ok(maximum) } fn mean_abs_err_with_policy( @@ -556,17 +601,19 @@ where A: Float, { let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); - let pairs = collect_policy_pairs(self, other, policy, context)?; - if pairs.is_empty() { + let (sum, count) = fold_policy_pairs( + self, + other, + policy, + context, + A::zero(), + |sum, first, second| sum + (first - second).abs(), + )?; + if count == 0 { return Err(StatisticalError::EmptyInput { context }); } - let sum = pairs - .iter() - .copied() - .fold(A::zero(), |sum, (a, b)| sum + (a - b).abs()) - .to_f64() - .expect("failed cast from type A to f64"); - Ok(sum / pairs.len() as f64) + let sum = sum.to_f64().expect("failed cast from type A to f64"); + Ok(sum / count as f64) } fn mean_sq_err_with_policy( @@ -577,17 +624,23 @@ where where A: Float, { - let sum = self - .sq_l2_dist_with_policy(other, policy)? - .to_f64() - .expect("failed cast from type A to f64"); - let pairs = collect_policy_pairs( + let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); + let (sum, count) = fold_policy_pairs( self, other, policy, - StatisticalErrorContext::new(StatisticalOperation::Deviation), + context, + A::zero(), + |sum, first, second| { + let difference = first - second; + sum + difference * difference + }, )?; - Ok(sum / pairs.len() as f64) + if count == 0 { + return Err(StatisticalError::EmptyInput { context }); + } + let sum = sum.to_f64().expect("failed cast from type A to f64"); + Ok(sum / count as f64) } fn root_mean_sq_err_with_policy( @@ -618,15 +671,18 @@ where private_impl! {} } -fn collect_policy_pairs( +fn fold_policy_pairs( first: &ArrayRef, second: &ArrayRef, policy: NumericPolicy, context: StatisticalErrorContext, -) -> Result, StatisticalError> + initial: R, + mut fold: F, +) -> Result<(R, usize), StatisticalError> where A: Float, D: Dimension, + F: FnMut(R, A, A) -> R, { if first.is_empty() { return Err(StatisticalError::EmptyInput { context }); @@ -639,17 +695,20 @@ where }); } - first + let mut state = initial; + let mut count = 0; + for (index, (first, second)) in first .iter() .copied() .zip(second.iter().copied()) .enumerate() - .try_fold(Vec::new(), |mut pairs, (index, (first, second))| { - if let Some(pair) = - apply_numeric_policy_pair(first, second, policy, context.with_index(index))? - { - pairs.push(pair); - } - Ok(pairs) - }) + { + if let Some((first, second)) = + apply_numeric_policy_pair(first, second, policy, context.with_index(index))? + { + state = fold(state, first, second); + count += 1; + } + } + Ok((state, count)) } diff --git a/src/entropy.rs b/src/entropy.rs index b12b967..183b373 100644 --- a/src/entropy.rs +++ b/src/entropy.rs @@ -3,7 +3,7 @@ use crate::errors::{ EmptyInput, MultiInputError, ShapeMismatch, StatisticalError, StatisticalErrorContext, StatisticalOperation, }; -use crate::policies::{apply_numeric_policy_pair, collect_with_numeric_policy, NumericPolicy}; +use crate::policies::{apply_numeric_policy_pair, fold_with_numeric_policy, NumericPolicy}; use ndarray::{Array, ArrayRef, Dimension, Zip}; use num_traits::Float; @@ -241,17 +241,23 @@ where A: Float, { let context = StatisticalErrorContext::new(StatisticalOperation::Entropy); - let values = collect_with_numeric_policy(self.iter().copied(), policy, context)?; - if values.is_empty() { + let (sum, count) = fold_with_numeric_policy( + self.iter().copied(), + policy, + context, + A::zero(), + |sum, value| { + if value == A::zero() { + sum + } else { + sum + value * value.ln() + } + }, + )?; + if count == 0 { return Err(StatisticalError::EmptyInput { context }); } - Ok(-values.iter().copied().fold(A::zero(), |sum, value| { - if value == A::zero() { - sum - } else { - sum + value * value.ln() - } - })) + Ok(-sum) } fn kl_divergence_with_policy( @@ -263,17 +269,17 @@ where A: Float, { let context = StatisticalErrorContext::new(StatisticalOperation::KLDivergence); - let pairs = collect_policy_pairs(self, q, policy, context)?; - if pairs.is_empty() { - return Err(StatisticalError::EmptyInput { context }); - } - Ok(-pairs.iter().copied().fold(A::zero(), |sum, (p, q)| { + let (sum, count) = fold_policy_pairs(self, q, policy, context, A::zero(), |sum, p, q| { if p == A::zero() { sum } else { sum + p * (q / p).ln() } - })) + })?; + if count == 0 { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(-sum) } fn cross_entropy_with_policy( @@ -285,31 +291,34 @@ where A: Float, { let context = StatisticalErrorContext::new(StatisticalOperation::CrossEntropy); - let pairs = collect_policy_pairs(self, q, policy, context)?; - if pairs.is_empty() { - return Err(StatisticalError::EmptyInput { context }); - } - Ok(-pairs.iter().copied().fold(A::zero(), |sum, (p, q)| { + let (sum, count) = fold_policy_pairs(self, q, policy, context, A::zero(), |sum, p, q| { if p == A::zero() { sum } else { sum + p * q.ln() } - })) + })?; + if count == 0 { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(-sum) } private_impl! {} } -fn collect_policy_pairs( +fn fold_policy_pairs( p: &ArrayRef, q: &ArrayRef, policy: NumericPolicy, context: StatisticalErrorContext, -) -> Result, StatisticalError> + initial: R, + mut fold: F, +) -> Result<(R, usize), StatisticalError> where A: Float, D: Dimension, + F: FnMut(R, A, A) -> R, { if p.shape() != q.shape() { return Err(StatisticalError::ShapeMismatch { @@ -319,17 +328,15 @@ where }); } - p.iter() - .copied() - .zip(q.iter().copied()) - .enumerate() - .try_fold(Vec::new(), |mut pairs, (index, (p, q))| { - if let Some(pair) = apply_numeric_policy_pair(p, q, policy, context.with_index(index))? - { - pairs.push(pair); - } - Ok(pairs) - }) + let mut state = initial; + let mut count = 0; + for (index, (p, q)) in p.iter().copied().zip(q.iter().copied()).enumerate() { + if let Some((p, q)) = apply_numeric_policy_pair(p, q, policy, context.with_index(index))? { + state = fold(state, p, q); + count += 1; + } + } + Ok((state, count)) } #[cfg(test)] diff --git a/src/policies.rs b/src/policies.rs index 247a449..3e8857d 100644 --- a/src/policies.rs +++ b/src/policies.rs @@ -166,13 +166,43 @@ where A: Float, I: IntoIterator, { - values - .into_iter() - .enumerate() - .try_fold(Vec::new(), |mut collected, (index, value)| { - if let Some(value) = apply_numeric_policy(value, policy, context.with_index(index))? { - collected.push(value); - } - Ok(collected) - }) + let values = values.into_iter(); + let (lower_bound, upper_bound) = values.size_hint(); + let mut collected = Vec::with_capacity(upper_bound.unwrap_or(lower_bound)); + + for (index, value) in values.enumerate() { + if let Some(value) = apply_numeric_policy(value, policy, context.with_index(index))? { + collected.push(value); + } + } + + Ok(collected) +} + +/// Folds values after applying the shared policy without materializing the +/// retained values. The returned count is the number of values included in the +/// calculation after omission. +pub(crate) fn fold_with_numeric_policy( + values: I, + policy: NumericPolicy, + context: StatisticalErrorContext, + initial: R, + mut fold: F, +) -> Result<(R, usize), StatisticalError> +where + A: Float, + I: IntoIterator, + F: FnMut(R, A) -> R, +{ + let mut state = initial; + let mut count = 0; + + for (index, value) in values.into_iter().enumerate() { + if let Some(value) = apply_numeric_policy(value, policy, context.with_index(index))? { + state = fold(state, value); + count += 1; + } + } + + Ok((state, count)) } diff --git a/src/summary_statistics/means.rs b/src/summary_statistics/means.rs index 03ebead..43bf1a8 100644 --- a/src/summary_statistics/means.rs +++ b/src/summary_statistics/means.rs @@ -4,7 +4,9 @@ use crate::errors::{ EmptyInput, MultiInputError, ShapeMismatch, StatisticalError, StatisticalErrorContext, StatisticalOperation, SummaryStatisticsError, }; -use crate::policies::{apply_numeric_policy_pair, collect_with_numeric_policy, NumericPolicy}; +use crate::policies::{ + apply_numeric_policy_pair, collect_with_numeric_policy, fold_with_numeric_policy, NumericPolicy, +}; use ndarray::{Array, ArrayBase, ArrayRef, Axis, Data, Dimension, Ix1, RemoveAxis}; use num_integer::IterBinomial; use num_traits::{Float, FromPrimitive, Zero}; @@ -105,8 +107,18 @@ where A: Float + FromPrimitive, { let context = StatisticalErrorContext::new(StatisticalOperation::Mean); - let values = collect_with_numeric_policy(self.iter().copied(), policy, context)?; - mean_values(&values, context) + let (sum, count) = fold_with_numeric_policy( + self.iter().copied(), + policy, + context, + A::zero(), + |sum, value| sum + value, + )?; + if count == 0 { + Err(StatisticalError::EmptyInput { context }) + } else { + Ok(sum / A::from_usize(count).expect("length must convert to A")) + } } fn weighted_mean_with_policy( @@ -118,8 +130,21 @@ where A: Float + FromPrimitive, { let context = StatisticalErrorContext::new(StatisticalOperation::WeightedMean); - let pairs = collect_weight_pairs(self, weights, policy, context)?; - weighted_mean_values(&pairs, context) + let ((weighted_sum, weight_sum), count) = fold_weighted_pairs( + self, + weights, + policy, + context, + (A::zero(), A::zero()), + |(weighted_sum, weight_sum), value, weight| { + (weighted_sum + value * weight, weight_sum + weight) + }, + )?; + if count == 0 { + Err(StatisticalError::EmptyInput { context }) + } else { + Ok(weighted_sum / weight_sum) + } } fn weighted_sum_with_policy( @@ -131,8 +156,19 @@ where A: Float + FromPrimitive, { let context = StatisticalErrorContext::new(StatisticalOperation::WeightedSum); - let pairs = collect_weight_pairs(self, weights, policy, context)?; - weighted_sum_values(&pairs, context) + let (sum, count) = fold_weighted_pairs( + self, + weights, + policy, + context, + A::zero(), + |sum, value, weight| sum + value * weight, + )?; + if count == 0 { + Err(StatisticalError::EmptyInput { context }) + } else { + Ok(sum) + } } fn weighted_mean_axis_with_policy( @@ -145,22 +181,32 @@ where A: Float + FromPrimitive, D: RemoveAxis, { - let context = - StatisticalErrorContext::axis(StatisticalOperation::WeightedMean, axis.index()); - let pairs = weighted_axis_values(self, axis, weights, policy, context)?; + let pairs = fold_weighted_axis( + self, + axis, + weights, + policy, + StatisticalOperation::WeightedMean, + (A::zero(), A::zero()), + |(weighted_sum, weight_sum), value, weight| { + (weighted_sum + value * weight, weight_sum + weight) + }, + )?; let shape = self.raw_dim().remove_axis(axis); let values = pairs .into_iter() .enumerate() - .map(|(lane, pairs)| { - weighted_mean_values( - &pairs, - StatisticalErrorContext::axis_lane( - StatisticalOperation::WeightedMean, - axis.index(), - lane, - ), - ) + .map(|(lane, ((weighted_sum, weight_sum), count))| { + let context = StatisticalErrorContext::axis_lane( + StatisticalOperation::WeightedMean, + axis.index(), + lane, + ); + if count == 0 { + Err(StatisticalError::EmptyInput { context }) + } else { + Ok(weighted_sum / weight_sum) + } }) .collect::, _>>()?; Ok(Array::from_shape_vec(shape, values) @@ -177,22 +223,30 @@ where A: Float + FromPrimitive, D: RemoveAxis, { - let context = - StatisticalErrorContext::axis(StatisticalOperation::WeightedSum, axis.index()); - let pairs = weighted_axis_values(self, axis, weights, policy, context)?; + let pairs = fold_weighted_axis( + self, + axis, + weights, + policy, + StatisticalOperation::WeightedSum, + A::zero(), + |sum, value, weight| sum + value * weight, + )?; let shape = self.raw_dim().remove_axis(axis); let values = pairs .into_iter() .enumerate() - .map(|(lane, pairs)| { - weighted_sum_values( - &pairs, - StatisticalErrorContext::axis_lane( - StatisticalOperation::WeightedSum, - axis.index(), - lane, - ), - ) + .map(|(lane, (sum, count))| { + let context = StatisticalErrorContext::axis_lane( + StatisticalOperation::WeightedSum, + axis.index(), + lane, + ); + if count == 0 { + Err(StatisticalError::EmptyInput { context }) + } else { + Ok(sum) + } }) .collect::, _>>()?; Ok(Array::from_shape_vec(shape, values) @@ -204,15 +258,17 @@ where A: Float + FromPrimitive, { let context = StatisticalErrorContext::new(StatisticalOperation::HarmonicMean); - let values = collect_with_numeric_policy(self.iter().copied(), policy, context)?; - if values.is_empty() { + let (sum, count) = fold_with_numeric_policy( + self.iter().copied(), + policy, + context, + A::zero(), + |sum, value| sum + value.recip(), + )?; + if count == 0 { return Err(StatisticalError::EmptyInput { context }); } - let sum = values - .iter() - .copied() - .fold(A::zero(), |sum, value| sum + value.recip()); - Ok((A::from_usize(values.len()).expect("length must convert to A") / sum).recip()) + Ok(A::from_usize(count).expect("length must convert to A") / sum) } fn geometric_mean_with_policy(&self, policy: NumericPolicy) -> Result @@ -220,15 +276,17 @@ where A: Float + FromPrimitive, { let context = StatisticalErrorContext::new(StatisticalOperation::GeometricMean); - let values = collect_with_numeric_policy(self.iter().copied(), policy, context)?; - if values.is_empty() { + let (sum, count) = fold_with_numeric_policy( + self.iter().copied(), + policy, + context, + A::zero(), + |sum, value| sum + value.ln(), + )?; + if count == 0 { return Err(StatisticalError::EmptyInput { context }); } - let sum = values - .iter() - .copied() - .fold(A::zero(), |sum, value| sum + value.ln()); - Ok((sum / A::from_usize(values.len()).expect("length must convert to A")).exp()) + Ok((sum / A::from_usize(count).expect("length must convert to A")).exp()) } fn raw_moment_with_policy( @@ -326,8 +384,15 @@ where A: Float + FromPrimitive, { let context = StatisticalErrorContext::new(StatisticalOperation::WeightedVariance); - let pairs = collect_weight_pairs(self, weights, policy, context)?; - weighted_variance_values(&pairs, ddof, context) + let (state, count) = fold_weighted_pairs( + self, + weights, + policy, + context, + (A::zero(), A::zero(), A::zero()), + update_weighted_variance, + )?; + weighted_variance_state(state.0, state.1, state.2, count, ddof, context) } fn weighted_std_with_policy( @@ -339,7 +404,16 @@ where where A: Float + FromPrimitive, { - Ok(self.weighted_var_with_policy(weights, ddof, policy)?.sqrt()) + let context = StatisticalErrorContext::new(StatisticalOperation::WeightedStandardDeviation); + let (state, count) = fold_weighted_pairs( + self, + weights, + policy, + context, + (A::zero(), A::zero(), A::zero()), + update_weighted_variance, + )?; + Ok(weighted_variance_state(state.0, state.1, state.2, count, ddof, context)?.sqrt()) } fn weighted_var_axis_with_policy( @@ -353,16 +427,25 @@ where A: Float + FromPrimitive, D: RemoveAxis, { - let context = - StatisticalErrorContext::axis(StatisticalOperation::WeightedVariance, axis.index()); - let pairs = weighted_axis_values(self, axis, weights, policy, context)?; + let pairs = fold_weighted_axis( + self, + axis, + weights, + policy, + StatisticalOperation::WeightedVariance, + (A::zero(), A::zero(), A::zero()), + update_weighted_variance, + )?; let shape = self.raw_dim().remove_axis(axis); let values = pairs .into_iter() .enumerate() - .map(|(lane, pairs)| { - weighted_variance_values( - &pairs, + .map(|(lane, (state, count))| { + weighted_variance_state( + state.0, + state.1, + state.2, + count, ddof, StatisticalErrorContext::axis_lane( StatisticalOperation::WeightedVariance, @@ -387,9 +470,37 @@ where A: Float + FromPrimitive, D: RemoveAxis, { - Ok(self - .weighted_var_axis_with_policy(axis, weights, ddof, policy)? - .mapv_into(|value| value.sqrt())) + let pairs = fold_weighted_axis( + self, + axis, + weights, + policy, + StatisticalOperation::WeightedStandardDeviation, + (A::zero(), A::zero(), A::zero()), + update_weighted_variance, + )?; + let shape = self.raw_dim().remove_axis(axis); + let values = pairs + .into_iter() + .enumerate() + .map(|(lane, (state, count))| { + weighted_variance_state( + state.0, + state.1, + state.2, + count, + ddof, + StatisticalErrorContext::axis_lane( + StatisticalOperation::WeightedStandardDeviation, + axis.index(), + lane, + ), + ) + .map(|value| value.sqrt()) + }) + .collect::, _>>()?; + Ok(Array::from_shape_vec(shape, values) + .expect("weighted standard-deviation lanes must match the output shape")) } fn kurtosis_with_policy(&self, policy: NumericPolicy) -> Result @@ -739,15 +850,18 @@ where Ok(sum / A::from_usize(values.len()).expect("length must convert to A")) } -fn collect_weight_pairs( +fn fold_weighted_pairs( values: &ArrayRef, weights: &ArrayRef, policy: NumericPolicy, context: StatisticalErrorContext, -) -> Result, StatisticalError> + initial: R, + mut fold: F, +) -> Result<(R, usize), StatisticalError> where A: Float, D: Dimension, + F: FnMut(R, A, A) -> R, { if values.shape() != weights.shape() { return Err(StatisticalError::ShapeMismatch { @@ -757,32 +871,40 @@ where }); } - values + let mut state = initial; + let mut count = 0; + for (index, (value, weight)) in values .iter() .copied() .zip(weights.iter().copied()) .enumerate() - .try_fold(Vec::new(), |mut pairs, (index, (value, weight))| { - if let Some(pair) = - apply_numeric_policy_pair(value, weight, policy, context.with_index(index))? - { - pairs.push(pair); - } - Ok(pairs) - }) + { + if let Some((value, weight)) = + apply_numeric_policy_pair(value, weight, policy, context.with_index(index))? + { + state = fold(state, value, weight); + count += 1; + } + } + Ok((state, count)) } -fn weighted_axis_values( +fn fold_weighted_axis( values: &ArrayRef, axis: Axis, weights: &ArrayRef, policy: NumericPolicy, - context: StatisticalErrorContext, -) -> Result>, StatisticalError> + operation: StatisticalOperation, + initial: R, + mut fold: F, +) -> Result, StatisticalError> where A: Float, D: Dimension, + R: Clone, + F: FnMut(R, A, A) -> R, { + let context = StatisticalErrorContext::axis(operation, axis.index()); if values.is_empty() { return Err(StatisticalError::EmptyInput { context }); } @@ -799,56 +921,56 @@ where .into_iter() .enumerate() .map(|(lane, values)| { - values + let context = StatisticalErrorContext::axis_lane(operation, axis.index(), lane); + let mut state = initial.clone(); + let mut count = 0; + for (index, (value, weight)) in values .iter() .copied() .zip(weights.iter().copied()) .enumerate() - .try_fold(Vec::new(), |mut pairs, (index, (value, weight))| { - if let Some(pair) = apply_numeric_policy_pair( - value, - weight, - policy, - StatisticalErrorContext::axis_lane(context.operation, axis.index(), lane) - .with_index(index), - )? { - pairs.push(pair); - } - Ok(pairs) - }) + { + if let Some((value, weight)) = + apply_numeric_policy_pair(value, weight, policy, context.with_index(index))? + { + state = fold(state, value, weight); + count += 1; + } + } + Ok((state, count)) }) .collect() } -fn weighted_sum_values( - pairs: &[(A, A)], - context: StatisticalErrorContext, -) -> Result +fn update_weighted_variance((weight_sum, mean, sum): (A, A, A), value: A, weight: A) -> (A, A, A) where A: Float, { - if pairs.is_empty() { - return Err(StatisticalError::EmptyInput { context }); - } - Ok(pairs - .iter() - .copied() - .fold(A::zero(), |sum, (value, weight)| sum + value * weight)) + let next_weight_sum = weight_sum + weight; + let value_minus_mean = value - mean; + let next_mean = mean + (weight / next_weight_sum) * value_minus_mean; + let next_sum = sum + weight * value_minus_mean * (value - next_mean); + (next_weight_sum, next_mean, next_sum) } -fn weighted_mean_values( - pairs: &[(A, A)], +fn weighted_variance_state( + weight_sum: A, + _mean: A, + sum: A, + count: usize, + ddof: A, context: StatisticalErrorContext, ) -> Result where A: Float, { - let weighted_sum = weighted_sum_values(pairs, context)?; - let weight_sum = pairs - .iter() - .copied() - .fold(A::zero(), |sum, (_, weight)| sum + weight); - Ok(weighted_sum / weight_sum) + if count == 0 { + return Err(StatisticalError::EmptyInput { context }); + } + if ddof.is_nan() || ddof < A::zero() || ddof > A::one() { + return Err(StatisticalError::InvalidDegreesOfFreedom { context }); + } + Ok(sum / (weight_sum - ddof)) } fn moments_values( @@ -916,35 +1038,6 @@ where } } -fn weighted_variance_values( - pairs: &[(A, A)], - ddof: A, - context: StatisticalErrorContext, -) -> Result -where - A: Float + FromPrimitive, -{ - if pairs.is_empty() { - return Err(StatisticalError::EmptyInput { context }); - } - let zero = A::zero(); - let one = A::one(); - if ddof < zero || ddof > one { - return Err(StatisticalError::InvalidDegreesOfFreedom { context }); - } - - let mut weight_sum = zero; - let mut mean = zero; - let mut sum = zero; - for &(value, weight) in pairs { - weight_sum = weight_sum + weight; - let value_minus_mean = value - mean; - mean = mean + (weight / weight_sum) * value_minus_mean; - sum = sum + weight * value_minus_mean * (value - mean); - } - Ok(sum / (weight_sum - ddof)) -} - /// Private function for `weighted_var` without conditions and asserts. fn inner_weighted_var( arr: &ArrayRef, diff --git a/tests/numeric_policy.rs b/tests/numeric_policy.rs index fba410d..02d5673 100644 --- a/tests/numeric_policy.rs +++ b/tests/numeric_policy.rs @@ -1,5 +1,5 @@ use approx::assert_abs_diff_eq; -use ndarray::{arr2, array}; +use ndarray::{arr2, array, Axis}; use ndarray_stats::{ errors::{StatisticalError, StatisticalOperation}, policies::NumericPolicy, @@ -92,6 +92,209 @@ fn infinity_remains_distinct_from_missing_values() { )); } +#[test] +fn policy_summary_methods_cover_scalar_and_axis_weighted_results() { + let data = array![1.0, f64::NAN, 3.0]; + let weights = array![1.0, f64::NAN, 1.0]; + let policy = NumericPolicy::omit_missing(); + + assert_abs_diff_eq!( + data.harmonic_mean_with_policy(policy).unwrap(), + 1.5, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + data.geometric_mean_with_policy(policy).unwrap(), + 3.0f64.sqrt(), + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + data.raw_moments_with_policy(2, policy).unwrap()[2], + 5.0, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + data.central_moments_with_policy(2, policy).unwrap()[2], + 1.0, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + data.standardized_moments_with_policy(2, policy).unwrap()[2], + 1.0, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + data.skewness_with_policy(policy).unwrap(), + 0.0, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + data.kurtosis_with_policy(policy).unwrap(), + 1.0, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + data.weighted_var_with_policy(&weights, 0.0, policy) + .unwrap(), + 1.0, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + data.weighted_std_with_policy(&weights, 0.0, policy) + .unwrap(), + 1.0, + epsilon = 1e-12 + ); + + let matrix = arr2(&[[1.0, f64::NAN, 3.0], [4.0, 5.0, 6.0]]); + let axis_weights = array![1.0, 2.0, 1.0]; + let means = matrix + .weighted_mean_axis_with_policy(Axis(1), &axis_weights, policy) + .unwrap(); + let sums = matrix + .weighted_sum_axis_with_policy(Axis(1), &axis_weights, policy) + .unwrap(); + let variances = matrix + .weighted_var_axis_with_policy(Axis(1), &axis_weights, 0.0, policy) + .unwrap(); + let standard_deviations = matrix + .weighted_std_axis_with_policy(Axis(1), &axis_weights, 0.0, policy) + .unwrap(); + + assert_abs_diff_eq!(means, array![2.0, 5.0], epsilon = 1e-12); + assert_abs_diff_eq!(sums, array![4.0, 20.0], epsilon = 1e-12); + assert_abs_diff_eq!(variances, array![1.0, 0.5], epsilon = 1e-12); + assert_abs_diff_eq!( + standard_deviations, + array![1.0, 0.5f64.sqrt()], + epsilon = 1e-12 + ); +} + +#[test] +fn policy_rank_and_entropy_methods_cover_omitted_observations() { + let data = arr2(&[[1.0, f64::NAN, 3.0], [2.0, 4.0, 6.0]]); + let policy = NumericPolicy::omit_missing(); + + let spearman = data.spearman_correlation_with_policy(policy).unwrap(); + let kendall = data.kendall_tau_with_policy(policy).unwrap(); + assert_abs_diff_eq!(spearman[[0, 1]], 1.0, epsilon = 1e-12); + assert_abs_diff_eq!(kendall[[0, 1]], 1.0, epsilon = 1e-12); + + let p = array![0.5, f64::NAN, 0.5]; + let q = array![0.5, 0.25, 0.25]; + let expected_cross_entropy = -(0.5 * 0.5f64.ln() + 0.5 * 0.25f64.ln()); + assert_abs_diff_eq!( + p.cross_entropy_with_policy(&q, policy).unwrap(), + expected_cross_entropy, + epsilon = 1e-12 + ); +} + +#[test] +fn policy_deviations_fold_without_changing_omission_semantics() { + let first = array![1.0, f64::NAN, 3.0]; + let second = array![1.0, 2.0, 2.0]; + let policy = NumericPolicy::omit_missing(); + + assert_eq!(first.count_eq_with_policy(&second, policy).unwrap(), 1); + assert_eq!(first.count_neq_with_policy(&second, policy).unwrap(), 1); + assert_abs_diff_eq!( + first.sq_l2_dist_with_policy(&second, policy).unwrap(), + 1.0, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + first.l2_dist_with_policy(&second, policy).unwrap(), + 1.0, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + first.l1_dist_with_policy(&second, policy).unwrap(), + 1.0, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + first.linf_dist_with_policy(&second, policy).unwrap(), + 1.0, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + first.mean_abs_err_with_policy(&second, policy).unwrap(), + 0.5, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + first.mean_sq_err_with_policy(&second, policy).unwrap(), + 0.5, + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + first.root_mean_sq_err_with_policy(&second, policy).unwrap(), + 0.5f64.sqrt(), + epsilon = 1e-12 + ); + assert_abs_diff_eq!( + first + .peak_signal_to_noise_ratio_with_policy(&second, 1.0, policy) + .unwrap(), + 10.0 * (2.0f64).log10(), + epsilon = 1e-12 + ); +} + +#[test] +fn policy_reports_empty_and_invalid_degrees_of_freedom_consistently() { + let all_missing = array![f64::NAN, f64::NAN]; + let other = array![1.0, 2.0]; + let empty_eq = all_missing + .count_eq_with_policy(&other, NumericPolicy::omit_missing()) + .unwrap_err(); + assert!(matches!( + empty_eq, + StatisticalError::EmptyInput { context } + if context.operation == StatisticalOperation::Deviation + )); + let empty_neq = all_missing + .count_neq_with_policy(&other, NumericPolicy::omit_missing()) + .unwrap_err(); + assert!(matches!( + empty_neq, + StatisticalError::EmptyInput { context } + if context.operation == StatisticalOperation::Deviation + )); + + let data = array![1.0, 2.0]; + let invalid_covariance = data + .into_shape_with_order((1, 2)) + .unwrap() + .cov_with_policy(f64::NAN, NumericPolicy::propagate()) + .unwrap_err(); + assert!(matches!( + invalid_covariance, + StatisticalError::InvalidDegreesOfFreedom { context } + if context.operation == StatisticalOperation::Covariance + )); + + let invalid_variance = array![1.0, 2.0] + .weighted_var_with_policy(&array![1.0, 1.0], f64::NAN, NumericPolicy::propagate()) + .unwrap_err(); + assert!(matches!( + invalid_variance, + StatisticalError::InvalidDegreesOfFreedom { context } + if context.operation == StatisticalOperation::WeightedVariance + )); + + let invalid_standard_deviation = array![1.0, 2.0] + .weighted_std_with_policy(&array![1.0, 1.0], f64::NAN, NumericPolicy::propagate()) + .unwrap_err(); + assert!(matches!( + invalid_standard_deviation, + StatisticalError::InvalidDegreesOfFreedom { context } + if context.operation == StatisticalOperation::WeightedStandardDeviation + )); +} + fn crate_context( operation: StatisticalOperation, axis: usize,