SciPy/stats.kendalltau

用 Fenwick 树统计 Kendall 相关系数的不一致对

SciPy 在 Kendall 秩相关计算中用 Fenwick 树维护已扫描样本的秩频数,避免逐对比较所有样本。

v1.15.2检索与计算核验 2026-10-12
相关概念

工程问题#

Kendall 相关系数需要统计两个变量排序不一致的样本对。直接枚举 n 个样本的所有配对需要二次级比较;同时,两个变量中的相等值不能被误计为排序相反。

实际实现#

kendalltau 先把 y 压缩成从 1 开始的稠密秩,再按 x 稳定排序。_kendall_dis 用 Fenwick 树保存已处理 x 组的 y 秩频数,对当前 y 查询不大于它的前缀计数,从已处理样本数中相减得到不一致对数。相同 x 的一组先全部查询,再统一插入,避免把组内平局当作不一致。

处理流程#

  1. 对 y 排序并压缩为正整数秩
  2. 稳定排序 x,保留同一 x 组中的次序
  3. 用 Fenwick 前缀和计算此前样本中 y 更大的数量
  4. 同一 x 组查询完毕后批量更新秩频数
  5. 结合两侧及共同平局计数,计算 Kendall 相关系数和检验结果

源码片段#

固定版本 · 原文摘录
scipy/stats/_stats_py.pyv1.15.2 · L5568–L5578
查看原文件
5568    size = x.size5569    perm = np.argsort(y)  # sort on y and convert y to dense ranks5570    x, y = x[perm], y[perm]5571    y = np.r_[True, y[1:] != y[:-1]].cumsum(dtype=np.intp)55725573    # stable sort on x and convert x to dense ranks5574    perm = np.argsort(x, kind='mergesort')5575    x, y = x[perm], y[perm]5576    x = np.r_[True, x[1:] != x[:-1]].cumsum(dtype=np.intp)55775578    dis = _kendall_dis(x, y)  # discordant pairs
kendalltau

将 y 转为从 1 开始的稠密秩,按 x 稳定排序后调用 _kendall_dis;数据准备保证后续 Fenwick 下标与同 x 分组有效。

scipy/stats/_stats.pyxv1.15.2 · L84–L101
查看原文件
84def _kendall_dis(intp_t[:] x, intp_t[:] y):85    cdef:86        intp_t sup = 1 + np.max(y)87        # Use of `>> 14` improves cache performance of the Fenwick tree (see gh-10108)88        intp_t[::1] arr = np.zeros(sup + ((sup - 1) >> 14), dtype=np.intp)89        intp_t i = 0, k = 0, size = x.size, idx90        int64_t dis = 09192    with nogil:93        while i < size:94            while k < size and x[i] == x[k]:95                dis += i96                idx = y[k]97                while idx != 0:98                    dis -= arr[idx + (idx >> 14)]99                    idx = idx & (idx - 1)100101                k += 1
_kendall_dis · prefix query

dis 先增加已处理样本数 i,再减去 y 前缀频数;idx & (idx - 1) 逐步清除最低位的 1,合并前缀区间。相同 x 的样本在这里先查询,不立即写树。

scipy/stats/_stats.pyxv1.15.2 · L103–L110
查看原文件
103            while i < k:104                idx = y[i]105                while idx < sup:106                    arr[idx + (idx >> 14)] += 1107                    idx += idx & -idx108                i += 1109110    return dis
_kendall_dis · point update

一组查询结束后插入各个 y 秩;idx += idx & -idx 更新覆盖该位置的父区间,避免组内相同 x 被计为不一致。

先展示公开接口的数据准备和调用点,再分开给出前缀查询与单点更新。复杂度针对秩计数内核;包含排序时通常为 O(n log n),不把此界限推广到所有 p 值计算方法。

收益与取舍#

工程分析

带来的收益

  • 当不同 y 秩数为 m 时,每次前缀查询和频数更新只访问 O(log m) 个数组位置
  • 分组查询再更新将 x 平局排除,前缀计数包含相等 y,从而也不把 y 平局计为不一致对

付出的代价

  • 需要排序、秩压缩及 O(m) 计数数组;并非无预处理的流式统计
  • 这是不一致对的计数内核,完整 Kendall 检验还需处理缺失值、平局修正及 p 值计算

可比较的选择

枚举全部样本对实现直观,适合小样本校验;样本增长时比较数为 O(n²)。

归并排序统计逆序对同样可做到 O(n log n),但需要在归并计数与平局处理之间明确边界,不能直接用普通逆序对计数替代全部 Kendall 逻辑。

版本与适用范围#

对应 SciPy v1.15.2 stats.kendalltau 的非缺失值主路径。输入由 Python 包装层排序并转为稠密秩后才进入 Cython 内核;nan_policy="omit" 等路径可能转交其他实现。树状数组使用带间隔的物理下标布局,>> 14 不改变 Fenwick 前缀和及单点更新的数学语义。

参考来源#

源码、官方文档与原始论文
01
scipy · kendalltau

scipy/stats/_stats_py.py · kendalltau

源码v1.15.2
02
scipy · _kendall_dis · prefix query

scipy/stats/_stats.pyx · _kendall_dis · prefix query

源码v1.15.2
03
scipy · _kendall_dis · point update

scipy/stats/_stats.pyx · _kendall_dis · point update

源码v1.15.2
19 个匹配案例