writing a function for columnwise computation of sample skewness in Rcpp I'm in trouble using sqrt()-function. I know that sqrt(x) works for NumericVector types (tested it in a separate file), however in my code (where I'm trying to pass doubles) it does not work.
Here's my code:
#include <Rcpp.h>
#include <cmath>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector colSkew(NumericMatrix x) {
int nc = x.ncol();
int nr = x.nrow();
NumericVector colS(nc);
for(int i = 0; i < nc; i++){
double cMean = mean(x( _ ,i));
double xSq = 0;
double cSt = 0;
for(int j = 0; j < nr; j++){
xSq += std::pow(x(j,i), 2.0);
cSt += std::pow(x(j,i) - cMean, 3.0);
}
double colMsq = nr * std::pow(cMean, 2.0);
double cTT = sqrt((xSq - colMsq)) / (nr - 1);
double colVar = cTT / (nr - 1);
double colNew = nr * std::pow(colVar, 3);
colS[i] = cSt / colNew;
}
return(colS);
}
I tried std::sqrt()
instead, giving me the "call to sqrt is ambiguous" error. For that one, neither Getting: error C2668: 'sqrt' : ambiguous call to overloaded function, nor http://dirk.eddelbuettel.com/code/rcpp.examples.html's
inline static double sqrt_double( double x ){ return ::sqrt( x ); }
helped. (The latter one has itself either the "no matching function" or the "ambiguous" problem). The code itself compiles, but does not work properly (I can call the function but results are not fine).
Reproducible example:
# First load the cpp function, however you want to;
Create example data:
A = matrix(rchisq(1000, 5), nrow = 100)
library(timeDate)
skew.1 = apply(A, 2, skewness)
skew.2 = colSkew(A)
# A custom R-function which should mimic the cpp function
colSkew.r = function(x){
nc = ncol(x)
nr = nrow(x)
colS = numeric(nc)
cMean = colMeans(x)
xSq = colSums(x^2)
cSt = 0
for(i in 1:nc){
cSt = sum((x[,i]-cMean[i])^3)
colMsq = nr * cMean[i]^2
cTT = sqrt((xSq[i] - colMsq) / (nr - 1))
colNew = nr * cTT^3
colS[i] = cSt / colNew
}
return(colS)
}
skew.3 = colSkew.r(A)
To be brief, there are few recurrent and common beginner errors we can avoid:
do not include unneeded headers (mostly harmless, but Rcpp
already brings math headers)
do not use a globally flattened namespace Rcpp
especially after you get compilation errors on symbol visibility
With that and the required minimal changes it builds and runs (and returns something non-sensical but I leave that for you to tackle :)
There are other ways to sort this out as std::
symbols and Rcpp
symbols generally coexist quite happily with distinct signatures.
#include <Rcpp.h>
// [[Rcpp::export]]
Rcpp::NumericVector colSkew(Rcpp::NumericMatrix x) {
int nc = x.ncol();
int nr = x.nrow();
Rcpp::NumericVector colS(nc);
for(int i = 0; i < nc; i++){
double cMean = Rcpp::mean(x( Rcpp::_ ,i));
double xSq = 0;
double cSt = 0;
for(int j = 0; j < nr; j++){
xSq += std::pow(x(j,i), 2.0);
cSt += std::pow(x(j,i) - cMean, 3.0);
}
double colMsq = nr * std::pow(cMean, 2.0);
double cTT = std::sqrt((xSq - colMsq)) / (nr - 1);
double colVar = cTT / (nr - 1);
double colNew = nr * std::pow(colVar, 3);
colS[i] = cSt / colNew;
}
return(colS);
}
/*** R
set.seed(42)
colSkew( matrix(rnorm(100), 100, 1) )
*/
> sourceCpp("answer.cpp")
> set.seed(42)
> colSkew( matrix(rnorm(100), 100, 1) )
[1] -448250877
>