feat: Missing data and infinity policies

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::<f64>();
+                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::<f64>();
+                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<f32> = 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();

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<A>` 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.