r - 是否有比 Base R 中的 expand.grid 更快的配对比较方法？
问题描述
同事，
我想估计从一组中随机选择的项目比从不同组中随机选择的项目得分更高的可能性。这就是优势概率，有时称为通用语言效应大小，例如参见：https ://rpsychologist.com/d3/cohend/ 。如果我们接受数据是正态分布的，这可以通过代数解决（McGraw 和 Wong (1992, Psychological Bulletin, 111)，但我知道我的数据不是正态分布的，因此这样的估计是不可靠的。
我的解决方案是简单的蛮力。使用样本数据，我们将一组中的每个观察值与另一组中的每个观察值进行比较，计算 a>b 的频率并将任何关系减半。
我的第一次尝试是For If Else循环，看起来像这样：
# Generate sample data
alpha <- as.integer(runif(2000, min = 0, max = 100))
beta <- as.integer(runif(2000, min = 0, max = 100))
nPairs <- length(alpha) * length(beta)
#initialise vars
bigger <- 0
equal <- 0.0
smaller <- 0
# Loop
for (i in alpha) {
for (j in beta) {
if (i > j) {bigger <- bigger + 1}
else
if (i < j) {smaller <- smaller + 1}
else
{equal <- equal + 0.5}
}
}
# record Probability a > b
PaGTb <- (bigger + equal) / nPairs
这可以完成工作，但速度非常慢，尤其是在Shiny应用程序中。
一位工作同事提醒我，R是基于向量的，并建议了该expand.grid函数。所以我的第二次尝试看起来像这样：
# generating sample data
alpha <- as.integer(runif(2000, min = 0, max = 100))
beta <- as.integer(runif(2000, min = 0, max = 100))
nPairs <- length(alpha) * length(beta)
# Each observation in alpha is paired with each observation in beta
c <- expand.grid(a = alpha, b = beta)
# Count frequency of a > b
aGTb <- length(which(c$a > c$b))
aEQb <- length(which(c$a == c$b))
aGTb <- aGTb + (0.5 * aEQb)
# record Probability a > b
PaGTb <- aGTb / nPairs
速度明显加快了！
但不是快一个量子步骤。通过模拟，我发现基于向量的方法对于许多配对比较（数百万）更快，但 For-If-Else 方法对于相对较少的比较（数千）更快。下表总结了我的 iMac 上 3334*3000=10,002,000 配对比较的 1000 次复制的结果。
time_ForIfElse time_Vector Ratio_Vector/ForIfElse
Min. :0.3514 Min. :0.2835 Min. :0.2713
1st Qu.:0.3648 1st Qu.:0.2959 1st Qu.:0.7818
Median :0.3785 Median :0.3102 Median :0.8115
Mean :0.4037 Mean :0.3322 Mean :0.8309
3rd Qu.:0.4419 3rd Qu.:0.3678 3rd Qu.:0.8668
Max. :1.4124 Max. :0.5896 Max. :1.4726
因此，对于我正在使用的数据类型，基于向量的方法比我原来的方法快大约 20%。但我仍然认为可能有更快的方法来处理这个问题。
社区有什么想法吗？
解决方案
这是一个完整的基础解决方案，它依赖于table()@chinsoon12 并受其启发。
f_int_AUC <- function(a, b){
A <- table(a)
A_Freq = as.vector(A)
A_alpha = as.integer(names(A))
B <- table(b)
B_Freq = as.vector(B)
B_beta = as.integer(names(B))
bigger = 0
equal = 0
for (i in seq_along(A_alpha)){
bigger = bigger + sum(A_Freq[i] * B_Freq[A_alpha[i] > B_beta])
equal = equal + 0.5 * sum(A_Freq[i] * B_Freq[A_alpha[i] == B_beta])
}
(bigger + equal) / (length(a) * length(b))
}
将它用作函数很重要——它在 4,000 行上比未编译的循环快大约 8 倍。
Unit: milliseconds
expr min lq mean median uq max neval
base_int_AUC_fx 1.0187 1.03865 1.146774 1.10930 1.13230 5.3215 100
base_int_AUC 8.3071 8.43380 9.309865 8.53855 8.88005 40.3327 100
该函数的分析表明该table()调用正在减慢该解决方案的速度。为了解决这个问题，tabulate()现在并入以显着提高速度。缺点tabulate()是它没有命名它的垃圾箱。因此，我们需要对 bin 进行归一化（h/t 为 @chinsoon12 以获得额外 20% 的改进）：
f_int_AUC2 <- function(a,b) {
# tabulate does not include 0, therefore we need to
# normalize to positive integers.
m <- min(c(a,b))
A_Freq = tabulate(a - min(m) + 1)
B_Freq = tabulate(b - min(m) + 1)
# calculate the outer product.
out_matrix <- outer(A_Freq, B_Freq)
bigger = sum(out_matrix[lower.tri(out_matrix)])
equal = 0.5 * sum(diag(out_matrix))
(bigger + equal) / length(a) / length(b)
}
需要注意的一件事是使用table()或tabulate()做出假设并且在浮点数上效果不佳。根据@Aaron 的建议，我查看了wilcox.text代码并根据问题进行了一些简化。注意我从函数中提取了代码——wilcox.test()函数是 437 行代码——这只是其中的 4 行，它提供了显着的速度提升。
f_wilcox_AUC <- function(a, b){
# r <- data.table::frank(c(a, b)) #is much faster
r <- rank(c(a, b))
n.x <- as.integer(length(a))
n.y <- as.integer(length(b))
# from wilcox.test STATISTIC; I am not a statistician and do not follow
# see reference at end or Google "r wilcox.test code"
{sum(r[seq_along(a)]) - n.x * (n.x + 1) / 2} / n.x / n.y
}
表现：
4,000 x 4,000
Unit: microseconds
expr min lq mean median uq max neval
bigstats_prive 719 748 792 797 817 884 10
chinsoon_DT 5910 6050 6280 6210 6350 7180 10
base_int_AUC 1060 1070 2660 1130 1290 16300 10
base_int_AUC2 237 250 1430 256 266 12000 10
wilcox_base 1050 1050 1830 1060 1070 8790 10
wilcox_dt 500 524 2940 530 603 16800 10
wilcox_test 11300 11400 11900 11500 11700 13400 10
40,000 x 40,000
Unit: milliseconds
expr min lq mean median uq max neval
bigstats_prive 4.03 4.07 5.13 4.11 4.27 66.40 100
chinsoon_DT 6.62 7.01 7.88 7.44 7.75 19.20 100
base_int_AUC 4.53 4.60 5.96 4.68 4.81 93.40 100
base_int_AUC2 1.03 1.07 1.14 1.08 1.12 3.70 100
wilcox_base 13.70 13.80 14.10 13.80 14.00 26.50 100
wilcox_dt 1.87 2.00 2.38 2.11 2.23 6.25 100
wilcox_test 115.00 118.00 121.00 121.00 124.00 135.00 100
400,000 x 400,000
Unit: milliseconds
expr min lq mean median uq max neval
bigstats_prive 50.3 54.0 55.0 54.7 56.8 58.6 10
chinsoon_DT 19.8 20.9 22.6 22.7 23.7 27.3 10
base_int_AUC 43.6 45.3 53.3 47.8 49.6 108.0 10
base_int_AUC2 10.0 10.4 11.8 10.7 13.8 14.8 10
wilcox_base 252.0 254.0 271.0 258.0 293.0 316.0 10
wilcox_dt 19.4 20.8 22.1 22.6 23.6 24.2 10
wilcox_test 1440.0 1460.0 1480.0 1480.0 1500.0 1520.0 10
总之，使用基本函数将花费 210 毫秒来进行 4,000 X 4,000 比较的 1,000 次复制：
# A tibble: 7 x 13
expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc total_time result memory time gc
<bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl> <int> <dbl> <bch:tm> <list> <list> <list> <list>
1 bigstats_prive 680.9us 692.8us 1425. 410.76KB 8.60 994 6 697.45ms <dbl [1]> <df[,3] [12 x 3]> <bch:tm> <tibble [1,000 x 3]>
2 chinsoon_DT 5.55ms 6.02ms 166. 539.92KB 4.78 972 28 5.85s <dbl [1]> <df[,3] [82 x 3]> <bch:tm> <tibble [1,000 x 3]>
3 base_int_AUC 981.8us 1.02ms 958. 750.44KB 10.7 989 11 1.03s <dbl [1]> <df[,3] [606 x 3]> <bch:tm> <tibble [1,000 x 3]>
4 base_int_AUC2 203.7us 209.3us 4743. 454.36KB 33.4 993 7 209.38ms <dbl [1]> <df[,3] [22 x 3]> <bch:tm> <tibble [1,000 x 3]>
5 wilcox_base 1.03ms 1.04ms 959. 359.91KB 4.82 995 5 1.04s <dbl [1]> <df[,3] [11 x 3]> <bch:tm> <tibble [1,000 x 3]>
6 wilcox_dt 410.4us 444.5us 2216. 189.09KB 6.67 997 3 450.01ms <dbl [1]> <df[,3] [9 x 3]> <bch:tm> <tibble [1,000 x 3]>
7 wilcox_test 11.35ms 11.44ms 85.6 1.16MB 1.66 981 19 11.45s <dbl [1]> <df[,3] [58 x 3]> <bch:tm> <tibble [1,000 x 3]>
编辑：outer()请参阅以前的方法以了解使用和的低效方法data.table::CJ()
参考资料： https ://github.com/SurajGupta/r-source/blob/master/src/library/stats/R/wilcox.test.R
推荐阅读
- asp.net - 如何将 RadioButtonList 项绑定到 ControlParameter
- c# - c#中的索引超出范围错误，“找不到第0列”
- python - youtube_dl python 包可以自动更新吗？
- arangodb - Sapper & Polka 认证
- python - Python PEP 8 和函数定义的垂直距离
- python - pip 使用 alpine 映像在 docker 中安装 pandas
- python - Seaborn 散点图在某些图上反转 y 轴？
- powershell - Powershell csv仅在特定列中查找带有反斜杠的值并重命名它们
- arrays - I 系列 RPGLE 在多个数组中搜索公共值
- python - 如何在按下 tkinter 时更改按钮状态