diff --git a/benches/summary_statistics.rs b/benches/summary_statistics.rs index b1b607b..c8786d3 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; @@ -82,6 +84,16 @@ fn descriptive_statistics(c: &mut Criterion) { }) }); + group.bench_with_input(format!("policy/{len}"), &data, |b, data| { + b.iter(|| { + let summary = black_box( + data.descriptive_statistics_with_policy(NumericPolicy::default()) + .unwrap(), + ); + black_box(score_f64(&summary)); + }) + }); + group.bench_with_input(format!("repeated/{len}"), &data, |b, data| { b.iter(|| { let result = ( @@ -119,6 +131,17 @@ fn descriptive_statistics_axis(c: &mut Criterion) { }) }); + 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()) + .unwrap(), + ); + let score = summaries.iter().map(score_f64).sum::(); + black_box(score); + }) + }); + group.bench_with_input(format!("repeated_axis0/{len}"), &data, |b, data| { b.iter(|| { black_box(score_repeated_axis(data, Axis(0), &axis_zero_weights)); @@ -133,6 +156,17 @@ fn descriptive_statistics_axis(c: &mut Criterion) { }) }); + 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()) + .unwrap(), + ); + let score = summaries.iter().map(score_f64).sum::(); + black_box(score); + }) + }); + group.bench_with_input(format!("repeated_axis1/{len}"), &data, |b, data| { b.iter(|| { black_box(score_repeated_axis(data, Axis(1), &axis_one_weights)); @@ -156,6 +190,17 @@ fn descriptive_statistics_f32(c: &mut Criterion) { 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( + data.descriptive_statistics_with_policy(NumericPolicy::default()) + .unwrap(), + ); + black_box(score_f32(&summary)); + }) + }); } group.finish(); diff --git a/docs/SUMMARY_STATISTICS.md b/docs/SUMMARY_STATISTICS.md new file mode 100644 index 0000000..22d75e1 --- /dev/null +++ b/docs/SUMMARY_STATISTICS.md @@ -0,0 +1,195 @@ +# 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. + +`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 + +- 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..f842715 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,160 @@ 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))); + } + + 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 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(_) => {} + None => { + include = false; + } + } + } + if include { + retained_observations.push(observation); + } + } + + let retained = retained_observations.len(); + if retained == 0 { + return Err(StatisticalError::EmptyInput { context }); + } + Ok(Array2::from_shape_fn( + (n_variables, retained), + |(variable, index)| data[[variable, retained_observations[index]]], + )) +} + +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.is_nan() || 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..ea1f2a8 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,278 @@ where Ok(psnr) } + fn count_eq_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let (count, included) = fold_policy_pairs( + self, + other, + policy, + StatisticalErrorContext::new(StatisticalOperation::Deviation), + 0usize, + |count, first, second| count + usize::from(first == second), + )?; + if included == 0 { + Err(StatisticalError::EmptyInput { + context: StatisticalErrorContext::new(StatisticalOperation::Deviation), + }) + } else { + Ok(count) + } + } + + fn count_neq_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let (count, included) = fold_policy_pairs( + self, + other, + policy, + StatisticalErrorContext::new(StatisticalOperation::Deviation), + 0usize, + |count, first, second| count + usize::from(first != second), + )?; + if included == 0 { + Err(StatisticalError::EmptyInput { + context: StatisticalErrorContext::new(StatisticalOperation::Deviation), + }) + } else { + Ok(count) + } + } + + fn sq_l2_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + 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) + } + + fn l2_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + 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( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); + 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(sum) + } + + fn linf_dist_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); + 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(maximum) + } + + fn mean_abs_err_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::Deviation); + 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 = sum.to_f64().expect("failed cast from type A to f64"); + Ok(sum / count as f64) + } + + fn mean_sq_err_with_policy( + &self, + other: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + 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 }); + } + 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( + &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 fold_policy_pairs( + first: &ArrayRef, + second: &ArrayRef, + policy: NumericPolicy, + context: StatisticalErrorContext, + 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 }); + } + if first.shape() != second.shape() { + return Err(StatisticalError::ShapeMismatch { + context, + first_shape: first.shape().to_vec(), + second_shape: second.shape().to_vec(), + }); + } + + let mut state = initial; + let mut count = 0; + for (index, (first, second)) in first + .iter() + .copied() + .zip(second.iter().copied()) + .enumerate() + { + 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 4ba9972..183b373 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, fold_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,109 @@ where Ok(cross_entropy) } + fn entropy_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::Entropy); + 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(-sum) + } + + fn kl_divergence_with_policy( + &self, + q: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::KLDivergence); + 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( + &self, + q: &ArrayRef, + policy: NumericPolicy, + ) -> Result + where + A: Float, + { + let context = StatisticalErrorContext::new(StatisticalOperation::CrossEntropy); + 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 fold_policy_pairs( + p: &ArrayRef, + q: &ArrayRef, + policy: NumericPolicy, + context: StatisticalErrorContext, + 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 { + context, + first_shape: p.shape().to_vec(), + second_shape: q.shape().to_vec(), + }); + } + + 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)] 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..3e8857d 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,102 @@ 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, +{ + 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/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..43bf1a8 100644 --- a/src/summary_statistics/means.rs +++ b/src/summary_statistics/means.rs @@ -1,6 +1,12 @@ 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, 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}; @@ -14,7 +20,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 +53,472 @@ 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 (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( + &self, + weights: &Self, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::WeightedMean); + 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( + &self, + weights: &Self, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::WeightedSum); + 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( + &self, + axis: Axis, + weights: &ArrayRef, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + D: RemoveAxis, + { + 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, ((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) + .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 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, (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) + .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 (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 }); + } + Ok(A::from_usize(count).expect("length must convert to A") / sum) + } + + fn geometric_mean_with_policy(&self, policy: NumericPolicy) -> Result + where + A: Float + FromPrimitive, + { + let context = StatisticalErrorContext::new(StatisticalOperation::GeometricMean); + 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 }); + } + Ok((sum / A::from_usize(count).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 (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( + &self, + weights: &Self, + ddof: A, + policy: NumericPolicy, + ) -> Result + where + A: Float + FromPrimitive, + { + 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( + &self, + axis: Axis, + weights: &ArrayRef, + ddof: A, + policy: NumericPolicy, + ) -> Result, StatisticalError> + where + A: Float + FromPrimitive, + D: RemoveAxis, + { + 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, (state, count))| { + weighted_variance_state( + state.0, + state.1, + state.2, + count, + 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, + { + 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 + 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 +836,208 @@ 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 fold_weighted_pairs( + values: &ArrayRef, + weights: &ArrayRef, + policy: NumericPolicy, + context: StatisticalErrorContext, + 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 { + context, + first_shape: values.shape().to_vec(), + second_shape: weights.shape().to_vec(), + }); + } + + let mut state = initial; + let mut count = 0; + for (index, (value, weight)) in values + .iter() + .copied() + .zip(weights.iter().copied()) + .enumerate() + { + 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 fold_weighted_axis( + values: &ArrayRef, + axis: Axis, + weights: &ArrayRef, + policy: NumericPolicy, + 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 }); + } + 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)| { + 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() + { + 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 update_weighted_variance((weight_sum, mean, sum): (A, A, A), value: A, weight: A) -> (A, A, A) +where + A: Float, +{ + 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_variance_state( + weight_sum: A, + _mean: A, + sum: A, + count: usize, + ddof: A, + context: StatisticalErrorContext, +) -> Result +where + A: Float, +{ + 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( + 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) + } + } +} + /// 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..02d5673 --- /dev/null +++ b/tests/numeric_policy.rs @@ -0,0 +1,306 @@ +use approx::assert_abs_diff_eq; +use ndarray::{arr2, array, Axis}; +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) + )); +} + +#[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, + 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), + }) ); }