我需要对一个结果变量进行采样，给定一个具有逐行结果概率的矩阵。
set.seed(1010) #reproducibility
#create a matrix of probabilities
#three possible outcomes, 10.000 cases
probabilities <- matrix(runif(10000*3),nrow=10000,ncol=3)
probabilities <- probabilities / Matrix::rowSums(probabilities)
我能想到的最快的方法是将apply((和sample((结合起来。
#row-wise sampling using these probabilities
classification <- apply(probabilities, 1, function(x) sample(1:3, 1, prob = x))
然而，在我所做的工作中，这是计算瓶颈。你知道如何加快这段代码的速度/如何更有效地进行采样吗？
谢谢！
RLave关于Rcpp可能是最佳选择的评论很到位(sample()也需要RcppArmadillo(；我使用以下C++代码创建了这样一个函数：
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadilloExtensions/sample.h>
using namespace Rcpp;
// [[Rcpp::export]]
IntegerVector sample_matrix(NumericMatrix x, IntegerVector choice_set) {
int n = x.nrow();
IntegerVector result(n);
for ( int i = 0; i < n; ++i ) {
result[i] = RcppArmadillo::sample(choice_set, 1, false, x(i, _))[0];
}
return result;
}
然后，我通过在我的R会话中提供了该功能
Rcpp::sourceCpp("sample_matrix.cpp")
现在，我们可以在R中针对您最初的方法进行测试，以及使用purrr::map()和lapply():的其他建议
set.seed(1010) #reproducibility
#create a matrix of probabilities
#three possible outcomes, 10.000 cases
probabilities <- matrix(runif(10000*3),nrow=10000,ncol=3)
probabilities <- probabilities / Matrix::rowSums(probabilities)
probabilities_list <- split(probabilities, seq(nrow(probabilities)))
library(purrr)
library(microbenchmark)
microbenchmark(
apply = apply(probabilities, 1, function(x) sample(1:3, 1, prob = x)),
map = map(probabilities_list, function(x) sample(1:3, 1, prob = x)),
lapply = lapply(probabilities_list, function(x) sample(1:3, 1, prob = x)),
rcpp = sample_matrix(probabilities, 1:3),
times = 100
)
Unit: milliseconds
expr min lq mean median uq max neval
apply 307.44702 321.30051 339.85403 342.36421 350.86090 434.56007 100
map 254.69721 265.10187 282.85592 286.21680 295.48886 363.95898 100
lapply 249.68224 259.70178 280.63066 279.87273 287.10062 691.21359 100
rcpp 12.16787 12.55429 13.47837 13.81601 14.25198 16.84859 100
cld
c
b
b
a
节省的时间相当可观。
如果您愿意将probabilities放在list中，purrr::map或lapply似乎会快一点：
probabilities <- matrix(runif(10000*3),nrow=10000,ncol=3)
probabilities <- probabilities / Matrix::rowSums(probabilities)
probabilities_list <- split(probabilities, seq(nrow(probabilities)))
library(purrr)
set.seed(1010)
classification_list <- map(probabilities_list, function(x) sample(1:3, 1, prob = x))
set.seed(1010)
classification_list <- lapply(probabilities_list, function(x) sample(1:3, 1, prob = x))
基准：
microbenchmark::microbenchmark(
apply = {classification = apply(probabilities, 1, function(x) sample(1:3, 1, prob = x))},
map = {classification = map(probabilities_list, function(x) sample(1:3, 1, prob = x))},
lapply = {classification = lapply(probabilities_list, function(x) sample(1:3, 1, prob = x))},
times = 100
)
# Unit: milliseconds
# expr min lq mean median uq max neval
# apply 39.92883 42.59249 48.39247 45.03080 47.86648 94.39828 100
# map 35.54077 37.13866 42.19719 39.95046 41.56323 66.05167 100
#lapply 34.54861 36.48664 42.69512 39.20139 52.31494 59.29200 100
100.000病例
# Unit: milliseconds
# expr min lq mean median uq max neval
# apply 457.5310 520.4926 572.5974 552.1674 611.5640 957.3997 100
# map 391.4751 457.7326 488.3286 482.1459 512.2054 899.1380 100
#lapply 386.2698 443.6732 491.9957 475.4160 507.3677 868.6725 100
您可以考虑
vapply和- 并行化：
parallel::parApply
使用probabilities矩阵：
set.seed(1010) #reproducibility
#create a matrix of probabilities
#three possible outcomes, 10.000 cases
probabilities <- matrix(runif(10000*3), nrow=10000,ncol=3)
probabilities <- probabilities / Matrix::rowSums(probabilities)
classification <- apply(probabilities, 1, function(x) sample(1:3, 1, prob = x))
vapply
通过为FUN.VALUE指定类，您可能能够使它变得快速。
classification2 <- vapply(split(probabilities, 1:nrow(probabilities)),
function(x) sample(1:3, 1, prob = x),
FUN.VALUE = integer(1), USE.NAMES = FALSE)
head(classification2)
#> [1] 1 3 3 1 2 3
平行包装
benchmarkme::get_cpu()
#> $vendor_id
#> [1] "GenuineIntel"
#>
#> $model_name
#> [1] "Intel(R) Core(TM) i5-4288U CPU @ 2.60GHz"
#>
#> $no_of_cores
#> [1] 4
在上述环境中，
cl <- parallel::makeCluster(4)
doParallel::registerDoParallel(cl, cores = 4)
parApply()可以做apply()做的事。
classification3 <- parallel::parApply(cl, probabilities, 1, function(x) sample(1:3, 1, prob = x))
head(classification3)
#> [1] 2 2 2 2 3 3
三者比较，包括apply()溶液、
microbenchmark::microbenchmark(
question = { # yours
apply(probabilities, 1, function(x) sample(1:3, 1, prob = x))
},
vapp = {
vapply(split(probabilities, 1:nrow(probabilities)), function(x) sample(1:3, 1, prob = x), FUN.VALUE = integer(1), USE.NAMES = FALSE)
},
parr = {
parallel::parApply(cl, probabilities, 1, function(x) sample(1:3, 1, prob = x))
}
)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> question 49.93853 58.39965 65.05360 62.98119 68.28044 182.03267 100
#> vapp 44.19828 54.84294 59.47109 58.56739 62.05269 146.14792 100
#> parr 43.33227 48.16840 53.26599 50.87995 54.17286 98.67692 100
parallel::stopCluster(cl)