How to compute rowSums in rcpp

Question

I'm converting an R function into Rcpp, where I have used the R function rowSums, which appears to not be a valid sugar expression in Rcpp. I found code for an Rcpp version of rowSums here. But I'm getting

error: use of undeclared identifier

when I use rowSumsC() in my main Rcpp function.

Is there a easy fix?

Edit: The code

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

Answer 1

Rcpp officially added rowSum support in 0.12.8. Therefore, there is no need to use rowSumsC function devised by Hadley in Advanced R.

Having said this, there are a few issues with the code.

---

Rcpp presently does not support Matrix to Vector or Matrix to Matrix computations. (Support for the later may be added per #583, though if needed one should consider using RcppArmadillo or RcppEigen). Therefore, the following line is problematic:

w = w / rowSums(w)

To address this, first compute the rowSums and then standardize the matrix using a traditional for loop. Note: Looping in C++ is very fast unlike R.

NumericVector summed_by_row = rowSums(w);

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

---

Next, C++ indices begin at 0 not 1. Therefore, the following for loop is problematic:

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

The fix:

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

---

Lastly, the parameters of the function can be reduced as some of the values are not relevant or are overridden.

The function declaration goes from:

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

To:

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

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

  .....

---

Putting all of the above feedback together, we get the desired function.

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)

Output

Before: 
0.0470993 0.125384
0.0285671 0.160996

rowSum: 0.172483 0.189563
After: 
0.273066 0.661436
0.165623 0.849300