r - 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;

}

它返回与 R 中的 dbinom() 等效的对数似然:

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

有人对更快地评估二项式对数似然有一些指导吗？可能是更快的代码或来自概率论的一些 hack。谢谢!

Answer 1

您的实现看起来不错。由于 R 的 dbinom() 已经在高效的 C 代码中实现，您可能不会显着改进它。我确实看到了一些可能会产生细微差别的事情(当你多次这样做时，这可能会有所帮助):

• 您可以使用 [ii] 而不是 (ii) 来避免边界检查，因为这听起来像是您不需要这样做的情况担心这一点(即，这不会是用户调用的函数，它只会在您的 C++ 代码中调用，在这种情况下，您的对象的设置方式可能不会成为问题)
• 您可以通过引用而不是通过值传递(参见，例如 here)

因此，我添加了您的函数的以下版本:

// [[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