企业官网建设流程全解析

CuPy 统计函数参考指南：顺序统计量、均值方差、相关性与直方图的 GPU 实现与用法

【免费下载链接】cupyNumPy & SciPy for GPU项目地址: https://gitcode.com/GitHub_Trending/cu/cupy

CuPy 是一套面向 GPU 的 NumPy/SciPy 兼容数值计算库，本文以其 API 参考文档 statistics.rst 为骨架，系统讲解 CuPy 统计子模块（cupy._statistics）提供的四大类函数：顺序统计量、均值与方差、相关性分析、直方图统计。读完本文，你将掌握每个函数的参数语义、与 NumPy 的兼容性边界、设备同步注意事项，以及它们背后的 CUDA Kernel 与 CUB 加速实现原理，能够直接在自己的 GPU 计算任务中正确选用并写出高性能统计代码。

一、统计模块在 CuPy 中的定位

CuPy 的公共 API 设计与 NumPy 保持一致。在 statistics.rst 中，全部统计函数被划分为四个分组，与 NumPy 官方文档的routines.statistics章节一一对应：

• Order statistics（顺序统计量）：ptp、percentile、quantile
• Averages and variances（均值与方差）：median、average、mean、std、var、nanmedian、nanmean、nanstd、nanvar
• Correlating（相关性）：corrcoef、correlate、cov
• Histograms（直方图）：histogram、histogram2d、histogramdd、bincount、digitize

从源码结构看，这些函数分散在 cupy/_statistics 目录的四个文件中：order.py（顺序统计量）、meanvar.py（均值方差）、correlation.py（相关性）、histogram.py（直方图），并在 cupy/init.py 中被统一导出到cupy命名空间顶层，因此用法与 NumPy 完全一致：cupy.mean(x)、cupy.histogram(x)等。此外，amin/amax/nanmin/nanmax也在该模块中实现并作为cupy.min/cupy.max导出。

需要特别留意的是，statistics.rst 中以#注释掉的条目（nanpercentile、nanquantile、histogram_bin_edges）表示当前版本尚未提供这些 API，使用前应通过hasattr(cupy, 'nanpercentile')等方式确认，避免导入错误。

二、顺序统计量：ptp、percentile 与 quantile

2.1 ptp：峰值到峰值的取值范围

ptp（peak-to-peak）返回数组沿指定轴的最大值与最小值之差：

import cupy as cp a = cp.array([[4, 9, 2], [10, 5, 13], [1, 14, 7]]) cp.ptp(a) # 13（整个数组的 max - min） cp.ptp(a, axis=0) # array([ 9, 9, 11]) cp.ptp(a, axis=1) # array([7, 8, 13])

参数axis缺省时作用于展平后的数组；out可指定输出数组；keepdims=True时被归约的轴保留为长度 1 的轴。从 order.py 的源码看，ptp直接委托给ndarray.ptp方法。注意：当某个归约切片内含有 NaN 时，对应的ptp结果也是 NaN；文档同时注明，若启用 cuTENSOR 加速器，含 NaN 的归约轴输出值可能被折叠（collapsed），这是加速路径与逐元素路径的已知差异。

2.2 percentile 与 quantile：分位数计算的完整参数

percentile计算 q-th 百分位数（q 取值范围 0～100），quantile计算 q-th 分位数（q 取值范围 0～1）。两者共享同一个底层实现_quantile_unchecked，区别仅在于百分位数先把q除以 100，再校验取值区间：

cp.percentile(a, 50) # 中位数 cp.percentile(a, [25, 50, 75]) # 四分位数，q 支持 tuple/list/ndarray cp.quantile(a, 0.5) cp.quantile(a, 0.5, axis=0) cp.quantile(a, 0.5, method='higher', keepdims=True)

核心参数如下：

method 插值方法的完整集合：从 order.py 的分发逻辑看，支持linear（默认）、lower、higher、midpoint、nearest、inverted_cdf，以及基于 Hyndman & Fan (1996) 连续插值法的hazen（H&F type 5）、weibull（type 6）、median_unbiased（type 8）、normal_unbiased（type 9）。后四者通过源码顶部的_QUANTILE_PARAMS字典映射为对应的(alpha, beta)参数对：

_QUANTILE_PARAMS = { 'hazen': (0.5, 0.5), 'weibull': (0, 0), 'median_unbiased': (1/3, 1/3), 'normal_unbiased': (3/8, 3/8), }

尚未实现的方法：averaged_inverted_cdf、closest_observation、interpolated_inverted_cdf（NumPy 1.22 新增）当前会直接抛出ValueError；测试文件 test_order.py 中也将这三个方法注释为TODO(takagi) Not implemented。传入未知方法名同样抛ValueError。

GPU 上的实现方式：_quantile_unchecked先把归约轴上的数据排序（ap.sort(axis=axis)），将归约轴搬到最后一维并展平，然后按插值方法把q换算成虚拟下标indices；当需要线性插值时，会现场编译一个名为cupy_percentile_weightnening的ElementwiseKernel（见 order.py），对排序后的数组按idx_below = floor(idx)、weight_above = idx - idx_below做加权插值。因此分位数计算在 GPU 上等价于"排序 + 插值核"，对大型数组有良好吞吐。

三、均值与方差：mean / std / var / average / median 及 NaN 安全变体

3.1 mean / std / var：直接的 ndarray 方法委托

mean、std、var在 meanvar.py 中直接委托给ndarray的同名方法，共享以下参数语义：

• axis：int、int 序列或None，默认对展平数组计算；
• dtype：指定计算与输出的数据类型（如用float32累积）；
• out：输出数组；
• keepdims：保留被归约轴；
• std/var额外支持ddof（delta degrees of freedom），默认为 0（除以 N）；设为 1 时除以N-1，得到无偏样本方差。

x = cp.arange(6, dtype=cp.float32).reshape(2, 3) cp.mean(x) # 2.5 cp.mean(x, axis=1) # array([1., 4.]) cp.var(x, ddof=1) cp.std(x, axis=0, keepdims=True)

3.2 average：支持权重的加权平均

average在 meanvar.py 中实现，返回沿轴的加权平均：

cp.average(cp.array([1, 2, 3, 4]), weights=cp.array([4, 3, 2, 1])) # 2.0 data = cp.arange(6).reshape((3, 2)) cp.average(data, axis=1, weights=cp.array([1./4, 3./4]))

关键行为：

• weights=None时退化为mean；
• 权重数组与a形状不同时，要求axis必须显式指定且权重为一维数组，其长度须等于该轴长度；
• 所有权重之和为 0 时抛出ZeroDivisionError；
• 整数或布尔输入会通过numpy.promote_types提升到至少float64的result_dtype；
• returned=True时返回(avg, scl)元组，其中scl是权重和；
• 注意：提供weights时该函数可能触发设备同步（源码中以cupy.any(scl == 0.0)注释# synchronize!标出），混合 CPU/GPU 流水线中需留意。

3.3 median 与 nanmedian

median委托给 Cython 层的_statistics._median（见 meanvar.py），支持axis（含多轴）、out、overwrite_input、keepdims。nanmedian有一个值得注意的分流：当输入 dtype 属于浮点或复数（'efdFD'）时走专门的_nanmedian实现，否则直接复用median（整数数组不含 NaN，无需特殊处理）。

cp.median(cp.array([[10, 7, 4], [3, 2, 1]])) # 3.5 cp.median(cp.array([[10, 7, 4], [3, 2, 1]]), axis=0) # array([6.5, 4.5, 2.5]) cp.nanmedian(cp.array([1.0, cp.nan, 3.0])) # 1.0

3.4 NaN 安全变体：nanmean / nanvar / nanstd / nanmin / nanmax

这一组函数忽略 NaN 值参与计算。NaN 安全归约的核心实现位于 Cython 文件 cupy/_core/_routines_statistics.pyx：例如_nanmean_func通过create_reduction_func生成归约核，_nanvar_core使用ReductionKernel，其归约逻辑在模板函数nanvar_impl中实现——对 NaN 元素贡献 0，对有效元素贡献(x - mean) * (x - mean)，最后除以max(_count - ddof, 0LL)（见 _routines_statistics.pyx 附近的实现），并提供了复数版本的专用内核_nanvar_core_complex64/complex128。

Python 层的分流规则（meanvar.py）：

• nanmean/nanvar/nanstd对整数与布尔 dtype（dtype.kind in 'biu'）直接退化为普通mean/var/std，避免无谓的 NaN 检查；
• nanmin/nanmax（order.py）先调用 Cython 层归约，再通过content.isnan(res).any()检测是否存在全 NaN 切片，若有则发出RuntimeWarning('All-NaN slice encountered')并返回 NaN。

cp.nanmean(cp.array([[1, cp.nan], [3, 4]])) # 2.666... cp.nanstd(cp.array([1.0, cp.nan, 2.0, 4.0]), ddof=1) cp.nanmin(cp.array([1.0, cp.nan, 3.0])) # 1.0

设备同步警告：nanmin/nanmax（以及带权重的average、histogram系列）文档明确标注 "This function may synchronize the device"，因为全 NaN 切片检测需要把结果取回主机判断。在追求极致吞吐的循环中应评估其影响。

四、相关性分析：corrcoef / cov / correlate

4.1 cov：协方差矩阵

cov在 correlation.py 中实现，计算协方差矩阵：

x = cp.array([[0, 2], [1, 1], [2, 0]]).T cp.cov(x) # 2x2 协方差矩阵 cp.cov(x, rowvar=False)

参数说明：

• y：额外的变量观测集合，与a按行拼接；
• rowvar=True（默认）表示每行是一个变量、每列是一次观测；False则转置解释；
• bias=False时按(N-1)归一化（无偏估计），bias=True时按N归一化；ddof非None时覆盖bias隐含的默认值（ddof=1无偏、ddof=0简单平均）；
• fweights（整数频率权重）与aweights（观测向量权重）均要求为cupy.ndarray，且长度须等于观测数；两者同时给出时按w *= aweights相乘；
• dtype缺省时结果至少为float64精度（通过numpy.promote_types与float64提升）。

实现上，cov先做均值中心化X -= X.mean(axis=1)[:, None]，再通过矩阵乘法X.dot(X_T.conj()) / fact得到协方差矩阵，其中fact由ddof/bias与权重共同决定；当自由度fact <= 0时发出RuntimeWarning并置 0。注意：corrcoef中的bias与ddof参数已废弃（传入会触发DeprecationWarning且不生效）。

4.2 corrcoef：皮尔逊相关系数

corrcoef在 correlation.py 中基于cov实现：先求协方差矩阵，取其对角线的实部开方得到标准差，再对矩阵逐项归一化out /= stddev[:, None]; out /= stddev[None, :]，最后把实部（以及复数的虚部）clip 到[-1, 1]区间：

cp.corrcoef(x) cp.corrcoef(x, y) # y 为额外变量集

4.3 correlate：一维互相关

correlate计算两个一维序列的离散互相关，mode支持valid（默认）、same、full。实现上先通过_choose_conv_method在直接卷积与 FFT 卷积之间做选择，然后分别调用_dot_convolve或_fft_convolve，核心等价于convolve(a, v[::-1])（correlation.py）。输入为空或非一维数组时会抛出ValueError。

cp.correlate(cp.array([1, 2, 3]), cp.array([0, 1, 0.5]), mode='full')

五、直方图：histogram 系列、bincount 与 digitize

5.1 histogram：一维直方图

histogram在 histogram.py 中实现，返回(hist, bin_edges)元组：

cp.histogram(cp.arange(10), bins=5) cp.histogram(x, bins=cp.array([0., 2., 4., 6., 8., 10.])) # 显式 bin 边界 cp.histogram(x, range=(0, 10), density=True) cp.histogram(x, bins=10, weights=w) # 加权直方图

参数语义与 NumPy 对齐：

• bins：整数（等宽分箱个数，须为正数）或一维数组（单调递增的 bin 边界）；字符串形式的 bin 算法（如'auto'、'fd'）当前抛出NotImplementedError（源码注释only integer and array bins are implemented）；
• range：(min, max)二元组，缺省为(x.min(), x.max())；边界外取值被忽略，且range同时影响自动分箱宽度；
• density=True：返回概率密度bin_count / sample_count / bin_volume；
• weights：与x同形状的权重数组，支持可安全转换到 float/complex 的 dtype（如Decimal等对象 dtype 不支持）；
• 复数输入抛出NotImplementedError，布尔输入会先告警并转为uint8。

GPU 实现亮点：直方图的核心是二分查找分箱 + 原子累加。源码中预编译了两个ElementwiseKernel：_histogram_kernel与_weighted_histogram_kernel（histogram.py），每个元素在 GPU 上通过 while 循环二分定位 bin，再用atomicAdd累加计数。当启用 CUB 加速器且数据量不超过0x7fffffff时，优先调用cub.cub_histogram（cupy/cuda/cub.pyx 提供绑定），并针对 NumPy 的"最后一 bin 右闭"语义做修正：整数 bin 上界 +1，浮点 bin 上界用cupy.nextafter外推一位；HIP 平台上计数临时用uint64。CUB 路径失败时自动回退到 CuPy 内核。该函数与histogramdd/histogram2d均在文档中标注 "may synchronize the device"。

5.2 histogram2d 与 histogramdd：多维直方图

• histogram2d(x, y, bins=10, range=None, weights=None, density=None)：二维直方图，返回(H, xedges, yedges)。bins可为单个整数（两个维度共用）、二元组（各维分箱数）或一维数组序列（各维边界）；bins为cupy.ndarray时表示两维共用同一组边界。
• histogramdd(sample, bins=10, range=None, weights=None, density=False)：D 维直方图，返回(H, edges)。sample为(N, D)数组（每行一个 D 维坐标点），或(X, Y, Z, ...)形式的坐标序列（内部用cupy.stack组装）；bins可为标量、长度为 D 的序列或各维边界数组序列。

histogramdd的实现流程（histogram.py）：逐维用linspace生成等宽边界或用给定数组作边界 → 用cupy.searchsorted(edges[i], sample[:, i], side='right')计算每个样本落入的 bin 编号（源码注明刻意避开cupy.digitize以规避 NumPy gh-11022 问题）→ 对恰好落在右边界上的样本编号 -1 修正 →cupy.ravel_multi_index展平为扁平索引 →cupy.bincount(xy, weights, minlength=...)统计 → reshape 后裁掉每个维度的离群 bin（首尾各一），density=True时再除以各维 bin 宽度与总样本数。

5.3 bincount：非负整数计数

bincount(x, weights=None, minlength=None)统计非负整数数组每个取值的出现次数，输出长度等于max(cupy.max(x) + 1, minlength)：

cp.bincount(cp.array([0, 1, 1, 3, 2, 1, 7])) # array([1, 3, 1, 1, 0, 0, 0, 1]) cp.bincount(cp.array([0, 1, 1, 3]), weights=cp.array([0.5, 0.25, 0.75, 1.0]))

约束：x必须是一维非负整数数组（浮点 dtype 抛TypeError，出现负数抛ValueError）；weights形状须与x相同；minlength须非负。计数内核_bincount_kernel同样使用atomicAdd，无权重时优先走 CUB 路径（HIP 平台除外），有权重时使用_bincount_with_weight_kernel并以float64累积。空输入返回numpy.intp类型的零数组。该函数会触发设备同步（需要取回x.max()确定输出长度）。

5.4 digitize：样本所属 bin 的下标

digitize(x, bins, right=False)返回x中每个值所属区间的下标，结果形状与x相同，bins须为一维单调数组：

x = cp.array([0.2, 6.4, 3.0, 1.6]) bins = cp.array([0.0, 1.0, 2.5, 4.0, 10.0]) cp.digitize(x, bins) # array([1, 4, 3, 2])

right=True时区间左开右闭（即bins[i-1] <= x < bins[i]），否则左闭右开（bins[i-1] < x <= bins[i]）。实现上直接转发给排序模块的_searchsorted（histogram.py）。与 NumPy 的一个刻意差异在文档注释中注明：为避免设备同步，digitize不会对bins的单调性做校验，非单调数组的结果是未定义的。

六、性能相关：加速器（CUB）与设备同步清单

统计模块的性能与行为受两个机制影响：

1. 例程加速器（routine accelerators）。histogram与bincount会遍历_accelerator.get_routine_accelerators()（由CUPY_ACCELERATORS环境变量或cupy._core.set_routine_accelerators配置），优先调用 CUB 的DeviceHistogram绑定 cub_histogram，失败或未启用时回退到 CuPy 自研atomicAdd内核。CUB 路径对输入规模有<= 0x7fffffff的限制（源码TODO(leofang): support >= 2^31 elements in x?）。

2. 设备同步点。以下函数在特定条件下会把设备数据取回主机或执行同步，密集调用时需权衡：

七、如何查阅与验证：文档、源码与测试的对应关系

如果你需要深入某个函数的边界行为，可以按"参考文档 → Python 封装 → Cython 内核 → 测试"的链路查阅：

1. API 参考：docs/source/reference/statistics.rst 列出全部函数及其分组，并链接到各函数生成的独立页面；
2. Python 封装层：cupy/_statistics 下的order.py、meanvar.py、correlation.py、histogram.py包含完整的 docstring、参数校验与 NumPy 兼容逻辑；
3. Cython 归约内核：cupy/_core/_routines_statistics.pyx 定义了cupy_min/cupy_max归约函数以及nanmean/nanvar/nanstd等 NaN 安全归约内核；
4. 测试套件：tests/cupy_tests/statistics_tests 下的test_order.py、test_meanvar.py、test_correlation.py、test_histogram.py覆盖了各函数的 dtype 组合、axis/keepdims/method组合与异常路径，是验证预期行为的第一手资料。例如 test_order.py 中列出了percentile/quantile全部受测的method取值，可作为"哪些插值方法可用"的权威清单。

八、使用建议与已知限制

• 与 NumPy 的双向互操作：以上函数均接受并返回cupy.ndarray，可与 NumPy 数组通过cupy.asarray/numpy.asarray互转；多数函数保持与 NumPy 相同的 dtype 提升规则（如cov/corrcoef至少float64精度）。
• 已知未实现项：nanpercentile、nanquantile、histogram_bin_edges（文档中已注释）；histogram的字符串 bin 算法；percentile/quantile的averaged_inverted_cdf、closest_observation、interpolated_inverted_cdf三种 method；cov的对象 dtype 权重与负数权重不报错的性能优化行为也与 NumPy 略有差异。
• 同步代价意识：在 GPU 流水线中尽量一次性提交大批量数据，避免在小数组上频繁调用会同步的统计函数，以免多次 kernel 启动与设备同步吃掉并行收益。

总而言之，CuPy 的统计模块在 API 层面完整对齐 NumPy，在实现层面则针对 GPU 特性做了排序 + 插值内核、原子累加直方图、CUB 加速回退与复数/NaN 专用归约内核等深度优化。将 statistics.rst 与上述源码、测试配合使用，你就能在 GPU 上写出既符合 NumPy 习惯又充分释放硬件性能的统计计算代码。

【免费下载链接】cupyNumPy & SciPy for GPU项目地址: https://gitcode.com/GitHub_Trending/cu/cupy