CuPy实现的两种置换测试方法结果差异原因排查求助

核心差异根源 #

你遇到的问题本质是两种"置换"方法生成的协变量矩阵分布完全不同，只有当测试统计量无全局偏移时，这种差异的影响被掩盖，一旦加入固定偏移（+10），差异会被急剧放大：

• 原方法cp.random.permutation(covariates)：本质是打乱协变量矩阵的行顺序，因此每个列的1的数量与原矩阵完全一致，严格保留了原协变量的列边缘分布。
• 替代方法cp.random.rand(*covariates.shape) < proportion_of_ones：每个元素独立生成，列中1的数量服从二项分布（期望等于原列和，但方差不为0），完全破坏了原协变量的列边缘分布。

偏移导致差异的具体原因 #

当测试统计量加入+10偏移时，单个协变量的统计量可拆解为：

sum(score * cov_col) = 10 * sum(cov_col) + sum(noise * cov_col)

• 原置换方法中，sum(cov_col)是固定值（与原列和一致），统计量的波动仅来自sum(noise * cov_col)；
• 替代方法中，sum(cov_col)是随机变量，其波动被10倍放大，直接导致统计量的整体分布与真实置换分布偏差极大，最终p值出现显著差异。

当去掉+10时，统计量仅为sum(noise * cov_col)，此时列和波动的影响被零均值的噪声弱化，两种方法的分布差异在大量置换后不足以导致p值出现显著偏差。

验证步骤 #

可以通过以下步骤确认上述结论：

1. 计算原协变量矩阵的列和：original_col_sums = cp.sum(covariates, axis=0)；
2. 生成一批原方法的置换矩阵，检查所有列和是否与original_col_sums完全一致；
3. 生成一批替代方法的矩阵，计算列和的均值、方差，对比与original_col_sums的差异；
4. 分别计算两种方法下统计量的分布（均值、方差、分位数），观察带偏移时的差异是否被放大。

高效且正确的解决方案 #

要兼顾性能与置换测试的正确性，需生成保留原列和的置换矩阵，推荐两种优化方案：

方案1：批量行置换（原方法的高效版） #

原方法的行置换本身是正确的，只需通过批量操作提升效率：

import cupy as cp

n_rows = 500000
n_cols = 2000
n_perms = 20000  # 批量生成置换数

# 预先生成批量随机行索引
perm_indices = cp.array([cp.random.permutation(n_rows) for _ in range(n_perms)])
# 批量置换行，得到(n_perms, n_rows, n_cols)的矩阵
permuted_covs = covariates[perm_indices]

# 批量计算统计量
score_with_shift = cp.random.randn(n_rows) + 10
perm_stats = cp.sum(score_with_shift[None, :, None] * permuted_covs, axis=1)

这种批量操作比单次循环置换效率提升数倍，且完全符合置换测试的要求。

方案2：按列比例生成（次优但高效） #

若必须使用独立生成的方式，需针对每个列使用对应列的1的比例，而非全局比例，能大幅降低差异：

# 计算每个列的1的比例
col_proportions = cp.sum(covariates, axis=0) / n_rows
# 按列生成布尔矩阵
rand_cov = cp.random.rand(n_rows, n_cols) < col_proportions[None, :]

注意：这种方法仍无法完全匹配置换分布（列内元素是独立的，而置换是依赖的），但比全局比例的结果更接近真实置换。

示例代码（可切换CuPy/Numpy） #

import cupy as cp  # 替换为import numpy as np即可在CPU运行

# 模拟数据
n_rows = 500000
n_cols = 2000
covariates = cp.random.rand(n_rows, n_cols) < 0.1
col_proportions = cp.sum(covariates, axis=0) / n_rows
global_proportion = cp.mean(covariates)

# 测试统计量
score_with_shift = cp.random.randn(n_rows) + 10
score_no_shift = cp.random.randn(n_rows)

# 参考统计量
ref_with_shift = cp.sum(score_with_shift[:, None] * covariates, axis=0)
ref_no_shift = cp.sum(score_no_shift[:, None] * covariates, axis=0)

# 批量行置换计算
n_perms = 1000
perm_indices = cp.array([cp.random.permutation(n_rows) for _ in range(n_perms)])
perm_covs = covariates[perm_indices]
stats_rows_with_shift = cp.sum(score_with_shift[None, :, None] * perm_covs, axis=1)
stats_rows_no_shift = cp.sum(score_no_shift[None, :, None] * perm_covs, axis=1)

# 全局比例生成计算
stats_rand_global_with_shift = []
stats_rand_global_no_shift = []
for _ in range(n_perms):
    rand_cov = cp.random.rand(n_rows, n_cols) < global_proportion
    stats_rand_global_with_shift.append(cp.sum(score_with_shift[:, None] * rand_cov, axis=0))
    stats_rand_global_no_shift.append(cp.sum(score_no_shift[:, None] * rand_cov, axis=0))
stats_rand_global_with_shift = cp.array(stats_rand_global_with_shift)
stats_rand_global_no_shift = cp.array(stats_rand_global_no_shift)

# 计算p值
p_rows_with_shift = (cp.sum(stats_rows_with_shift >= ref_with_shift, axis=0) + 1) / (n_perms + 1)
p_rand_global_with_shift = (cp.sum(stats_rand_global_with_shift >= ref_with_shift, axis=0) + 1) / (n_perms + 1)
p_rows_no_shift = (cp.sum(stats_rows_no_shift >= ref_no_shift, axis=0) + 1) / (n_perms + 1)
p_rand_global_no_shift = (cp.sum(stats_rand_global_no_shift >= ref_no_shift, axis=0) + 1) / (n_perms + 1)

# 输出差异
print("带偏移时p值差异均值:", cp.mean(cp.abs(p_rows_with_shift - p_rand_global_with_shift)))
print("无偏移时p值差异均值:", cp.mean(cp.abs(p_rows_no_shift - p_rand_global_no_shift)))

内容的提问来源于stack exchange，提问作者Dan Bolfter