Medians

作者：Libor Spacek

寻找中位数的新算法，100% 安全的 Rust 实现。

use medians::{*,algos::*};

简介

在统计学和数据分析中，寻找中位数是一个常见任务。至少它应该是这样的，因为中位数是比平均值更稳定的集中趋势度量。同样，mad（绝对差的中位数）是比标准差更稳定的数据分散度量，而标准差受平方异常值的影响。中位数和 mad 并没有被足够地使用，这主要是由于实践历史上的原因：它们更难计算。这里提出的快速算法为这种情况提供了一种补救措施。

我们在 rstats 中论证，使用几何中位数是多维数据特征化的最稳定方法。本仓库解决了这个一维情况。

请参阅 tests.rs 以获取使用示例。它们的自动生成输出也可以通过点击此文档顶部的“测试”图标，然后检查最新日志来找到。

使用概述

根据最终数据类型（即输入向量/切片中项的类型）选择最佳方法和函数。

• u8 -> 函数 medianu8
• u64 -> 函数 medianu64
• f64 -> trait Medianf64 的方法
• T 自定义可量化为 u64 -> trait Median 的 uqmedian 方法
• T 自定义可比较 -> trait Median 的 qmedian_by 方法
• T 自定义可比较但不量化 -> trait Median 的通用方法 median_by。

算法分析

短原始类型最好通过基数搜索处理。我们已经为 u8 和 u64

/// Medians of u8 end type by fast radix search
pub fn medianu8(s: &[u8]) -> Result<ConstMedians<u8>, Me>;
/// Medians of u64 end type by fast recursive radix search
pub fn medu64(s: &mut [u64]) -> Result<(u64, u64), Me>;

更复杂的数据类型需要通用比较搜索，请参阅 median_by。可以通过对数据列表进行排序然后选择其中点来简单地找到中位数。最佳比较排序算法的复杂度为 O(n*log(n))。然而，可能的更快的 median 算法具有 O(n) 的复杂度。它们基于观察，即数据需要全部完全排序，只需要分区和计数。因此，简单的排序方法无法竞争，自版本 2.0.0 起已被删除。

Floyd-Rivest (1975)：目前认为“中位数的中位数”是比较算法的“艺术水平”。它将数据分成五个项目的组，通过排序找到每个组的中位数，然后找到五个这些中位数的中位数，依此类推，直到只剩下一个。然后将其用作原始数据分区的枢轴。这样的枢轴会产生良好的分区，虽然不是完美的二等分。因此，计数和迭代仍然是必要的。

找到最佳可能的枢轴估计并不是主要目标。真正的目标是尽快（计数）消除（计数）异常数据项。因此，必须考虑到估计枢轴所花费的时间。可以满足不那么优化的枢轴，同时平均更快地找到中位数。无论如何，有效的分区是必须的。

设我们平均每次分区后剩余的项目比率为 rs，Floyd-Rivest 的为 rf。通常，1/2 <= rf <= rs < 1，即 rf 更优，更接近完美的二等分（比率为 1/2）。假设我们可以在 Floyd-Rivest 完成一个所需的时间内执行两个分区（因为它们的枢轴选择较慢）。那么，只要 rs^2 < rf，这完全可能，并且在实践中似乎也是如此。例如，rf=0.65（几乎最优），rs=0.8（非常次优），但 rs^2 < rf。尽管如此，投入一些与数据长度成比例的计算努力来选择枢轴是值得的。

我们引入了另一种新算法，作为函数 medianu64

/// Fast medians of u64 end type by binary partitioning
pub fn medianu64(s: &mut [u64]) -> Result<ConstMedians<u64>, Me>

在 u64 数据上实现，它比 median_by 的一般目的枢轴选择快约两倍。数据通过单个位值进行分区，完全避免了枢轴估计的开销。该算法通常收敛良好。然而，当数据恰好全部聚集在很小的值范围内时，它会变慢。

我们通用中值算法的主要功能总结

• 线性复杂度。
• 快速（就地）迭代分区为三个子范围（较小、等于、较大），最小化数据移动和内存管理。
• 简单的枢轴选择策略：三个样本的中位数（只需要三次比较）。在迭代过程中，真正糟糕的枢轴只偶尔出现。对于更长的数据，我们使用三个中位数的中位数。

特质 Medianf64

/// Fast 1D medians of floating point data, plus related methods
pub trait Medianf64 {
    /// Median of f64s, NaNs removed
    fn medf_checked(self) -> Result<f64, Me>;
    /// Median of f64s, including NaNs
    fn medf_unchecked(self) -> f64;
    /// Iterative weighted median
    fn medf_weighted(self, ws: Self, eps: f64) -> Result<f64, Me>;
    /// Zero mean/median data produced by subtracting the centre
    fn medf_zeroed(self, centre: f64) -> Vec<f64>;
    /// Median correlation = cosine of an angle between two zero median vecs
    fn medf_correlation(self, v: Self) -> Result<f64, Me>;
    /// Median of absolute differences (MAD).
    fn madf(self, centre: f64) -> f64;
}

特质 Median

这些方法专门用于通用、任意复杂和/或大型数据类型。在分区等过程中，数据永远不会被复制。

其大部分方法都接受一个比较闭包 c，该闭包返回其泛型参数之间的排序，泛型参数为 &T。这允许以任何数量的不同方式比较任何自定义类型。

其大部分方法都接受一个量化闭包 q，它将泛型参数转换为 f64。这不仅简化了标准 Rust as 和 .into() 转换，而且还提供了多种灵活的方式来量化更复杂的自定义数据类型。

使用较弱的局部序数比较代替数值比较。搜索算法保持不变。唯一的额外成本是额外的引用层，以防止数据复制。

median_by()

对于所有可以量化为 f64 的端类型，我们简单地平均偶数长度数据的两个中点来获得单个中位数（类型为 f64）。当数据项无法量化为 f64 时，这不再可能。然后应使用 median_by。它返回 Medians 枚举类型中的两个中间值，首先是较小的那个。

/// Enum for results of odd/even medians
pub enum Medians<'a, T> {
    /// Odd sized data results in a single median
    Odd(&'a T),
    /// Even sized data results in a pair of (centered) medians
    Even((&'a T, &'a T)),
}

/// Fast 1D generic medians, plus related methods
pub trait Median<'a, T> {
    /// Median by comparison `c`, at the end quantified to a single f64 by `q`
    fn qmedian_by(
        self,
        c: &mut impl FnMut(&T, &T) -> Ordering,
        q: impl Fn(&T) -> f64,
    ) -> Result<f64, Me>;
    /// Median of types quantifiable to u64 by `q`, at the end converted to a single f64.  
    /// For data that is already `u64`, use function `medianu64`
    fn uqmedian(
            self,
            q: impl Fn(&T) -> u64,
        ) -> Result<f64, Me>;
    /// Median by comparison `c`, returns odd/even result
    fn median_by(self, c: &mut impl FnMut(&T, &T) -> Ordering) -> Result<Medians<'a, T>, Me>;
    /// Zero mean/median data, produced by subtracting the centre
    fn zeroed(self, centre: f64, quantify: impl Fn(&T) -> f64) -> Result<Vec<f64>, Me>;
    /// Median correlation = cosine of an angle between two zero median Vecs
    fn med_correlation(
        self,
        v: Self,
        c: &mut impl FnMut(&T, &T) -> Ordering,
        q: impl Fn(&T) -> f64,
    ) -> Result<f64, Me>;
    /// Median of absolute differences (MAD).
    fn mad(self, centre: f64, quantify: impl Fn(&T) -> f64) -> f64;
}

发行说明

版本 3.0.12 - 添加更快的 medu64，偶数版本仍在进行中。修复了一个错误。

版本 3.0.11 - 向 Median 特质添加了方法 uqmedian，用于将类型通过某些闭包 q 量化为 u64。修复了 oddmedian_by 中的最近错误，其中枢轴引用未及时保存。

版本 3.0.10 - 添加了 medianu64。它在 u64 数据上比通用 median_by 更快。它使用一种新的按位分区算法，从而避免了枢轴估计的复杂性。

版本 3.0.9 - 改进了大数据集的枢轴估计。

版本 3.0.8 - 添加了 implementation.rs 模块并重新组织了源代码。

版本 3.0.7 - 添加了 medf_weighted，应用 &[f64] 权重。

版本 3.0.6 - 将 part、ref_vec 和 deref_vec 移动到 crate Indxvec 中，以允许它们更广泛地使用。

版本 3.0.5 - 废弃代码修剪。

版本 3.0.4 - 一些代码简化。

版本 3.0.3 - 更新 dev 依赖 ran 到 2.0。

版本 3.0.2 - 添加了函数 medianu8，它通过超级快速基数搜索找到中值字节。将添加更多原始类型。

版本 3.0.1 - 将 correlation 重命名为 med_correlation 以避免与其他地方的名称冲突。

版本 3.0.0 - 在速度、通用性和重命名方面进行了多项改进。