如果你真的关心效率,这段简短的Rcpp代码将难以被超越。将以下内容存储在文件中,比如说
/tmp/rnormClamp.cpp
:
#include <Rcpp.h>
using namespace Rcpp;
NumericVector rnormClamp(int N, int mi, int ma) {
NumericVector X = rnorm(N, 0, 1);
return clamp(mi, X, ma);
}
使用sourceCpp()
(同时也可以使用Rcpp)进行构建和运行。在我的电脑上,实际的绘制和夹紧大约需要4毫秒:
R> sourceCpp("/tmp/rnormClamp.cpp")
R> system.time(X <- rnormClamp(50000, -3, 3))
user system elapsed
0.004 0.000 0.004
R> summary(X)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.00000 -0.67300 -0.00528 0.00122 0.68500 3.00000
R>
clamp()
糖函数在 Romain 的 SO 回答中被介绍,并且指出您需要 Rcpp 的版本为 0.10.2。
编辑:根据 Ben 的提示,我似乎误解了。这是 C++ 和 R 的混合内容:
List rnormSelect(int N, int mi, int ma) {
RNGScope scope;
int N2 = N * 1.25;
NumericVector X = rnorm(N2, 0, 1);
LogicalVector ind = (X < mi) | (X > ma);
return List::create(X, ind);
}
可以将其附加到早期的文件中。然后:
R> system.time({ Z <- rnormSelect(50000, -3, 3);
+ X <- Z[[1]][ ! Z[[2]] ]; X <- X[1:50000]})
user system elapsed
0.008 0.000 0.009
R> summary(X)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.00000 -0.68200 -0.00066 -0.00276 0.66800 3.00000
R>
我将重新审视逻辑索引和行子集,这部分内容我需要查阅。也许明天吧。但9毫秒还算不错 :)
编辑2: 看起来我们确实没有逻辑索引。我们需要添加它。此版本是“手动”执行的,但比从R索引略快。
NumericVector rnormSelect2(int N, int mi, int ma) {
RNGScope scope;
int N2 = N * 1.25;
NumericVector X = rnorm(N2, 0, 1);
LogicalVector ind = (X >= mi) & (X <= ma);
NumericVector Y(N);
int k=0;
for (int i=0; i<N2 & k<N; i++) {
if (ind[i]) Y(k++) = X(i);
}
return Y;
}
输出结果:
R> system.time(X <- rnormSelect2(50000, -3, 3))
user system elapsed
0.004 0.000 0.007
R> summary(X)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.99000 -0.66900 -0.00258 0.00223 0.66700 2.99000
R> length(X)
[1] 50000
R>
truncnorm::rtruncnorm()
感兴趣。 - hplieninger