2013-02-14 2 views
2
)

R 코드의 길이를 L (psi) 및 행렬 (L, L)의 벡터를 사용하고 요소 별 연산을 수행하는 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)) 
    } 

CPP :

#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; 
} 
+0

R 코드에'apply'를 사용하지 말고'log' 변수와 함께'colSums'를 사용하여 제품을 얻을 수 있습니다. – James

+1

첫 번째'apply'는't (phi) * psi'와 같습니다. 더 빠를 것입니다. – James

+1

로그를 생각한 다음 exp'ing하고 summing하는 것이 prods보다 훨씬 빠릅니까? – scottyaz

답변

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보다 훨씬 빠릅니다. 출력이 모두 1의 두 번째 부분에서 1보다 작은 수의 곱셈의 만기 숫자이기 때문에 내가 사용한 테스트 데이터가 관련 있는지 확실하지 않습니다. 이러한 작은 숫자의 세부 사항을 기록하려면 로그 축척으로 작업하는 것이 좋습니다.