用 Fenwick 树统计 Kendall 相关系数的不一致对
SciPy 在 Kendall 秩相关计算中用 Fenwick 树维护已扫描样本的秩频数,避免逐对比较所有样本。
工程问题#
Kendall 相关系数需要统计两个变量排序不一致的样本对。直接枚举 n 个样本的所有配对需要二次级比较;同时,两个变量中的相等值不能被误计为排序相反。
实际实现#
kendalltau 先把 y 压缩成从 1 开始的稠密秩,再按 x 稳定排序。_kendall_dis 用 Fenwick 树保存已处理 x 组的 y 秩频数,对当前 y 查询不大于它的前缀计数,从已处理样本数中相减得到不一致对数。相同 x 的一组先全部查询,再统一插入,避免把组内平局当作不一致。
处理流程#
- 对 y 排序并压缩为正整数秩
- 稳定排序 x,保留同一 x 组中的次序
- 用 Fenwick 前缀和计算此前样本中 y 更大的数量
- 同一 x 组查询完毕后批量更新秩频数
- 结合两侧及共同平局计数,计算 Kendall 相关系数和检验结果
源码片段#
固定版本 · 原文摘录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将 y 转为从 1 开始的稠密秩,按 x 稳定排序后调用 _kendall_dis;数据准备保证后续 Fenwick 下标与同 x 分组有效。
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 += 1dis 先增加已处理样本数 i,再减去 y 前缀频数;idx & (idx - 1) 逐步清除最低位的 1,合并前缀区间。相同 x 的样本在这里先查询,不立即写树。
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一组查询结束后插入各个 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 前缀和及单点更新的数学语义。