如何在 rcpp 中计算 rowSums

How to compute rowSums in rcpp

我正在将 R 函数转换为 Rcpp,其中我使用了 R 函数 rowSums,它在 Rcpp 中似乎不是有效的糖表达式。我找到了 rowSums here 的 Rcpp 版本的代码。但是我得到

error: use of undeclared identifier

当我在我的主要 Rcpp 函数中使用 rowSumsC() 时。

有简单的解决方法吗?

编辑:代码

cppFunction(
  "NumericMatrix Expcpp(NumericVector x, NumericMatrix w,
   NumericVector mu, NumericVector var, NumericVector prob, int k) {
   for (int i=1; i<k; ++i){
   w(_,i) = prob[i] * dnorm(x,mu[i], sqrt(var[i]));
   }
   w = w / rowSums(w)
   return w;
}")

Rcpp officially added rowSum support in 0.12.8。因此,无需使用 Hadley 在 Advanced R 中设计的 rowSumsC 函数。

话虽如此,代码还是有一些问题。


Rcpp 目前 支持MatrixVectorMatrixMatrix 计算。 (可以根据 #583, though if needed one should consider using RcppArmadillo or RcppEigen 添加对后者的支持)。因此,下面一行是有问题的:

w = w / rowSums(w)

为了解决这个问题,首先计算 rowSums,然后使用传统的 for 循环对矩阵进行标准化。 注意: C++ 中的循环非常快,不像 R

NumericVector summed_by_row = rowSums(w);

for (int i = 0; i < k; ++i) {
  w(_,i) = w(_,i) / summed_by_row[i];
}

接下来,C++ 索引从 0 而不是 1 开始。因此,下面的for循环是有问题的:

for (int i=1; i<k; ++i)

修复:

for (int i=0; i<k; ++i)

最后,可以减少函数的参数,因为某些值不相关或被覆盖。

函数声明来自:

NumericMatrix Expcpp(NumericVector x, NumericMatrix w,
   NumericVector mu, NumericVector var, NumericVector prob, int k)

收件人:

NumericMatrix Expcpp(NumericVector x, NumericVector mu, NumericVector var, NumericVector prob) {

  int n = x.size();
  int k = mu.size();
  NumericMatrix w = no_init(n,k); 

  .....

综合以上所有反馈,我们得到了想要的功能。

Rcpp::cppFunction(
  'NumericMatrix Expcpp(NumericVector x, NumericVector mu, NumericVector var, NumericVector prob) {

  int n = x.size();
  int k = mu.size();

  NumericMatrix w = no_init(n,k); 

  for (int i = 0; i < k; ++i) { // C++ indices start at 0
     w(_,i) = prob[i] * dnorm(x, mu[i], sqrt(var[i]));
  }

  Rcpp::Rcout << "Before: " << std::endl << w << std::endl;

  NumericVector summed_by_row = rowSums(w);

  Rcpp::Rcout << "rowSum: " << summed_by_row << std::endl;

  // normalize by column to mimic R
  for (int i = 0; i < k; ++i) {
    w(_,i) = w(_,i) / summed_by_row[i];
  }

  Rcpp::Rcout << "After: " << std::endl << w << std::endl;

  return w;
  }')

set.seed(51231)
# Test values
n <- 2
x <- seq_len(n)
mu <- x
var <- x
prob <- runif(n)

mat <- Expcpp(x, mu, var, prob)

输出

Before: 
0.0470993 0.125384
0.0285671 0.160996

rowSum: 0.172483 0.189563
After: 
0.273066 0.661436
0.165623 0.849300