2

我正在尝试使用长度为 L (psi) 的向量和维度矩阵 (L,L) 并执行一些元素操作的代码来加速一些R代码。Rcpp有没有更有效的方法来使用 Rcpp 进行这些元素操作？

r：

 UpdateLambda <- function(psi,phi){
                                      # updated full-day infection probabilites
    psi.times.phi <- apply(phi,1,function(x) x*psi)
    ## return Lambda_{i,j} = 1 - \prod_{j} (1 - \psi_{i,j,t} \phi_{i,j})
    apply(psi.times.phi,2,function(x) 1-prod(1-x))
  }

cp:

#include <Rcpp.h>
#include <algorithm>
using namespace Rcpp;


// [[Rcpp::export]]
NumericVector UpdateLambdaC(NumericVector psi,
                NumericMatrix phi
                ){

  int n = psi.size();
  NumericMatrix psi_times_phi(n,n);
  NumericVector tmp(n,1.0);
  NumericVector lambda(n);

  for(int i=0; i<n;i++){
    psi_times_phi(i,_) = psi*phi(i,_);
   }

  for(int i=0; i<n;i++){
     // \pi_{j} (1- \lambda_{i,j,t})
    for(int j=0; j<n;j++){
      tmp[i] *= 1-psi_times_phi(i,j);
    }
     lambda[i] = 1-tmp[i];
  }

  return lambda;
}
4

1 回答 1

1

您可以apply用矢量化替代品替换您的循环。

第一个相当于：

t(phi)*psi

第二个：

1-exp(colSums(log(1-psi.times.phi)))

#test data
phi <- matrix(runif(1e6),1e3)
psi <- runif(1e3)

#new function
UpdateLambda2 <- function(psi,phi) 1-exp(colSums(log(1-t(phi)*psi)))

#sanity check
identical(UpdateLambda(psi,phi),UpdateLambda2(psi,phi))
[1] TRUE

#timings
library(rbenchmark)
benchmark(UpdateLambda(psi,phi),UpdateLambda2(psi,phi))
                     test replications elapsed relative user.self sys.self
1  UpdateLambda(psi, phi)          100   16.05    1.041     15.06     0.93
2 UpdateLambda2(psi, phi)          100   15.42    1.000     14.19     1.19

好吧，看起来它并没有太大的区别，这非常令人惊讶，因为colSums它通常比apply. 我不确定我使用的测试数据是否相关，因为输出是 all1的第二部分中小于 1 的乘法数。如果您想注意这些小数字的细节，无论如何您最好使用对数刻度。

于 2013-02-14T14:11:36.523 回答