Rcpp中二项式似然的快速评估

Question

我需要快速评估大量的二项式似然。因此，我想在Rcpp中实现这个功能。其中一种方法是下面的。

#include <RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector eval_likelihood(arma::vec Yi,
                              arma::vec Ni,
                              arma::vec prob){

  // length of vector
  int N = prob.n_rows;

  // storage for evaluated log likelihoods
  NumericVector eval(N);

  for(int ii = 0; ii < N; ii++){

  int y = Yi(ii); // no. of successes
  int n = Ni(ii); // no. of trials
  double p = prob(ii); // success probability

  eval(ii) = R::dbinom(y,n,p,true); // argument 4 is set to true to return log-likelihood

  }

  return eval;

}

返回等价的对数似然 dbinom() 在R中。

Rcpp::sourceCpp("dbinom.cpp") #source Rcpp script

# fake data
Yi    = 1:999  
Ni    = 2:1000
probs = runif(999)

evalR    = dbinom(Yi, Ni, probs, log = T) # vectorized solution in R
evalRcpp = eval_likelihood(Yi, Ni, probs) # my Rcpp solution

identical(evalR,evalRcpp)
[1] TRUE

总的来说，这是一个不错的结果。但是，矢量化的R方案平均比我的天真Rcpp方案略快。

microbenchmark::microbenchmark(R    = dbinom(Yi, Ni, probs, log = T),
                               Rcpp = eval_likelihood(Yi, Ni, probs))

Unit: microseconds
 expr     min      lq     mean   median       uq      max neval cld
    R 181.753 182.181 188.7497 182.6090 189.4515  286.100   100   a
 Rcpp 178.760 179.615 197.5721 179.8285 184.7470 1397.144   100   a

有谁能指导我们如何更快的评估二项式对数似然？可能是更快的代码，也可能是概率论的一些黑客。谢谢！我需要评估大量的二项式对数似然。

Answer 1

你的实现看起来不错。作为R的 dbinom() 已经在高效的C代码中实现，你可能不会 大大 在此基础上进行改进。我确实看到了几件可能会产生微小差异的事情（当你做了很多次时，这可能会有帮助）。

• 你可以使用 [ii] 而非 (ii) 以避免边界检查，因为听起来你的情况是不需要担心这个问题的（也就是说，这不会是一个用户调用的函数，它只会在你的C++代码中被调用，而你的对象大概是以这样一种方式设置的，这不会是一个问题）。
• 你可以通过引用来传递，而不是通过值来传递（例如，见 这里的)

所以，我把你的函数添加到下面的版本中。

// [[Rcpp::export]]
NumericVector eval_likelihood2(const arma::vec& Yi,
                               const arma::vec& Ni,
                               const arma::vec& prob){

    // length of vector
    int N = prob.n_rows;

    // storage for evaluated log likelihoods
    NumericVector eval(N);

    for(int ii = 0; ii < N; ii++){

        int y = Yi[ii]; // no. of successes
        int n = Ni[ii]; // no. of trials
        double p = prob[ii]; // success probability

        eval[ii] = R::dbinom(y,n,p,1); // argument 4 is set to true to return log-likelihood

    }

    return eval;

}

你可以看到，我只是改变了这两样东西。

我也用了稍微大一点的数据作为基准，不过我也为你原来的小例子添加了基准。

Rcpp::sourceCpp("so.cpp") #source Rcpp script

# fake data
Yi    = 1:99999
Ni    = 2:100000
probs = runif(99999)

evalR     = dbinom(Yi, Ni, probs, log = T) # vectorized solution in R
evalRcpp  = eval_likelihood(Yi, Ni, probs) # my Rcpp solution
evalRcpp2 = eval_likelihood(Yi, Ni, probs) # my Rcpp solution

identical(evalR,evalRcpp)
# [1] TRUE
identical(evalR,evalRcpp2)
# [1] TRUE

microbenchmark::microbenchmark(R     = dbinom(Yi, Ni, probs, log = T),
                               Rcpp  = eval_likelihood(Yi, Ni, probs),
                               Rcpp2 = eval_likelihood2(Yi, Ni, probs))

Unit: milliseconds
  expr      min       lq     mean   median       uq      max neval
     R 7.427669 7.577011 8.565015 7.650762 7.916891 62.63154   100
  Rcpp 7.368547 7.858408 8.884823 8.014881 8.353808 63.48417   100
 Rcpp2 6.952519 7.256376 7.859609 7.376959 7.829000 12.51065   100

Yi    = 1:999
Ni    = 2:1000
probs = runif(999)
microbenchmark::microbenchmark(R     = dbinom(Yi, Ni, probs, log = T),
                               Rcpp  = eval_likelihood(Yi, Ni, probs),
                               Rcpp2 = eval_likelihood2(Yi, Ni, probs))

Unit: microseconds
  expr    min       lq     mean   median       uq     max neval
     R 90.073 100.5035 113.5084 109.5230 122.5260 188.304   100
  Rcpp 90.188  97.8565 112.9082 105.2505 122.4255 172.975   100
 Rcpp2 86.093  92.0745 103.9474  97.9380 113.2660 148.591   100