Умножение мудрых матричных элементов в Rcpp

Я пытаюсь ускорить код R с Rcpp, который принимает вектор длины L (psi) и матрицу измерений (L, L) и выполняет некоторые элементарные операции. Есть ли более эффективный способ выполнять эти элементарные операции с Rcpp?Умножение мудрых матричных элементов в Rcpp

г:

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)) 
    }

каст:

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

источник

2013-02-14 scottyaz

Похоже, вы можете избежать использования 'apply' в своем R-коде и использовать' colSums' с 'log'ed-переменными для получения продуктов. – James

Первый 'apply' эквивалентен' t (phi) * psi', который должен быть быстрее – James

, вы думаете, что журнал, а затем exp'ing и suming будут намного быстрее, чем prods? – scottyaz

Вы можете заменить вам apply петли с vectorised альтернатив.

Первый эквивалентно:

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 во второй части. В любом случае, вы можете лучше работать в масштабе журнала, если хотите отметить детали таких небольших чисел.

источник

2013-02-14 14:11:36 James

Умножение мудрых матричных элементов в Rcpp

ответ

Смежные вопросы