2012-10-23 7 views
12

R에서 데이터 행렬이 100 x 10 행렬 X이고 가능한 값 (0, 1, 2, 3)을 갖는 100 요소 벡터 t가있는 경우, 우리가 쉽게 간단한 구문을 사용하여 X의 행렬 y를 찾을 수 있습니다 문제는,논리 문장과 일치하는 Rcpp 행렬의 서브 세트

y = X[t == 1, ] 

을하지만, 내가 어떻게 할 수 Rcpp의 NumericMatrix와?
(또는 더 일반적으로, 내가 할 수있는 방법이 C++의 모든 컨테이너를?에서) 더크의 힌트에

덕분에, 그것은

NumericMatrix X(dataX); 
IntegerVector T(dataT); 
mat Xmat(X.begin(), X.nrow(), X.ncol(), false); 
vec tIdx(T.begin(), T.size(), false); 
mat y = X.rows(find(tIdx == 1)); 

내가 원하는 일을 할 수있는 것 같다,하지만 너무 긴 것 같다.

답변

7

내가 아는 가장 가까운 RcppArmadillo을 통해 액세스 할 수 Armadillosubmat() 기능과 결합 find() 함수의 조합입니다.

편집 : 이것은 패치를 통해 추가 할 수있는 과정입니다. 아무도이 문제를 해결할 동기가 충분하지 않으면 rcpp-devel 메일 링리스트로 오십시오.

+0

예,이를 추가하면 개발 및 테스트 꽤 많이 필요합니다. 따라서 전용 자금이 제공되지 않는 한 빨리 발생하지는 않습니다. –

6

설탕으로보고 싶습니다. 불행히도, 나는 그것을 구현할 수있는 자격이 없다. 여기에는 여전히 여러 가지 다른 솔루션이 있습니다.

첫째, 나는 (이 작업을 얻을 수 colvec 대신 vectIdxXmat.rows(... 대신 X.rows(... 징 이순신 리아 코드를 일부 수정했습니다 :

mat Xmat(X.begin(), X.nrow(), X.ncol(), false); 
colvec tIdx(T.begin(), T.size(), false); 
mat y = Xmat.rows(find(tIdx == 1)); 

둘째, 여기에 세 가지 기능이 있습니다 논리 문을 기반으로 모든 부분 집합 행렬을 벤치 마크합니다. 함수는 arma 또는 rcpp 인수와 반환 값을가집니다. 두 개는 Gong-Yi Liao의 솔루션을 기반으로하고 하나는 간단한 루프 기반 솔루션을 사용합니다.

N (행) = 100, P (T == 1) = 0.3

   expr min  lq median  uq max 
1 submat_arma(X, T) 5.009 5.3955 5.8250 6.2250 28.320 
2 submat_arma2(X, T) 4.859 5.2995 5.6895 6.1685 45.122 
3 submat_rcpp(X, T) 5.831 6.3690 6.7465 7.3825 20.876 
4  X[T == 1, ] 3.411 3.9380 4.1475 4.5345 27.981 

N (행) = 10000, P (T == 1) = 0.3

   expr  min  lq median  uq  max 
1 submat_arma(X, T) 107.070 113.4000 125.5455 141.3700 1468.539 
2 submat_arma2(X, T) 76.179 80.4295 88.2890 100.7525 1153.810 
3 submat_rcpp(X, T) 244.242 247.3120 276.6385 309.2710 1934.126 
4  X[T == 1, ] 229.884 236.1445 263.5240 289.2370 1876.980 

submat.cpp는

#include <RcppArmadillo.h> 
// [[Rcpp::depends(RcppArmadillo)]] 

using namespace Rcpp; 
using namespace arma; 

// arma in; arma out 
// [[Rcpp::export]] 
mat submat_arma(arma::mat X, arma::colvec T) { 
    mat y = X.rows(find(T == 1)); 
    return y; 
} 

// rcpp in; arma out 
// [[Rcpp::export]] 
mat submat_arma2(NumericMatrix X, NumericVector T) { 
    mat Xmat(X.begin(), X.nrow(), X.ncol(), false); 
    colvec tIdx(T.begin(), T.size(), false); 
    mat y = Xmat.rows(find(tIdx == 1)); 
    return y; 
} 

// rcpp in; rcpp out 
// [[Rcpp::export]] 
NumericMatrix submat_rcpp(NumericMatrix X, LogicalVector condition) { 
    int n=X.nrow(), k=X.ncol(); 
    NumericMatrix out(sum(condition),k); 
    for (int i = 0, j = 0; i < n; i++) { 
     if(condition[i]) { 
      out(j,_) = X(i,_); 
      j = j+1; 
     } 
    } 
    return(out); 
} 


/*** R 
library("microbenchmark") 

# simulate data 
n=100 
p=0.3 
T=rbinom(n,1,p) 
X=as.matrix(cbind(rnorm(n),rnorm(n))) 

# compare output 
identical(X[T==1,],submat_arma(X,T)) 
identical(X[T==1,],submat_arma2(X,T)) 
identical(X[T==1,],submat_rcpp(X,T)) 

# benchmark 
microbenchmark(X[T==1,],submat_arma(X,T),submat_arma2(X,T),submat_rcpp(X,T),times=500) 

# increase n 
n=10000 
p=0.3 
T=rbinom(n,1,p) 
X=as.matrix(cbind(rnorm(n),rnorm(n))) 
# benchmark 
microbenchmark(X[T==1,],submat_arma(X,T),submat_arma2(X,T),submat_rcpp(X,T),times=500) 

*/