我有一个对称的概率矩阵,对角线条目为空。假设类似
0 0.5 0.1 0.6
0.5 0 0.2 0.1
0.1 0.2 0 0.2
0.6 0.1 0.2 0
我想绘制一个虚拟矩阵,使条目 [i,j] 的概率成为概率矩阵中条目 [i,j] 的概率。请注意,我的概率矩阵是一个 Armadillo 矩阵(一个大矩阵 5000x5000)。当然,对角线虚拟变量应该为空,因为它们的概率为空。我构建了两个函数来做到这一点,但它们并不快。我应该在循环中多次采样这个矩阵。
mat binom1(mat& prob){
int n=prob.n_rows;
mat sample(n,n,fill::zeros);
NumericVector temp(2);
for(int i(0);i<n-1;++i){
for(int j(i+1);j<n;++j){
temp=rbinom(2,1,prob(i,j));
sample(i,j)=temp(0); sample(j,i)=temp(1);
}
}
return sample;
}
mat binom2(mat& prob){
int n=prob.n_rows;
mat sample(n,n);
for(int i(0);i<n;++i){
for(int j(0);j<n;++j){
sample(i,j)=as<double>(rbinom(1,1,prob(i,j)));
}
}
return sample;
}
两者都比 R 中的矢量化 rbinom 慢。
z=matrix(runif(1000^2),1000) #just an example for 1000x1000 matrix
microbenchmark(rbinom(nrow(z)^2,1,z),binom1(z),binom2(z))
结果
expr min lq mean median uq max
rbinom(nrow(z)^2, 1, z) 95.43756 95.94606 98.29283 97.5273 100.3040 108.2293
binom1(z) 131.33937 133.25487 139.75683 136.4530 139.5511 229.0484
binom2(z) 168.38226 172.60000 177.95935 175.6447 180.9531 277.3501
有没有办法让代码更快?
我看到一个例子 here 。但就我而言,概率位于 Armadillo 矩阵中
最佳答案
鉴于几乎重复的答案,您可以使用:
mat binom3(const mat& prob) {
int n = prob.n_rows;
mat sample(n, n);
std::transform(prob.begin(), prob.end(), sample.begin(),
[=](double p){ return R::rbinom(1, p); });
return sample;
}
微基准测试:
Unit: milliseconds
expr min lq mean median uq max neval
rbinom(length(z), 1, z) 46.88264 47.28971 48.09543 47.66346 48.40734 65.29790 100
binom1(z) 76.98416 82.60813 84.93669 83.51432 84.04780 126.46992 100
binom2(z) 96.20707 98.59145 101.99215 99.56175 102.02750 153.04754 100
binom3(z) 34.01417 34.49066 35.12199 34.93946 35.47979 38.22539 100
关于r - 具有 Armadillo 矩阵概率的矢量化 Rcpp rbinom,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/46642128/