write a rolling sd function in R and implement in C-CodePudding

I have to write implement R function for C . for e.g I am trying to calculate rolling SD but below code is not working. any help will be highly appreciated. #Below code is working fine

library(roll)

n <- 150
x <- rnorm(n)
x
weights <- 0.9 ^ (n:1)
weights
roll_sd(x, width = 5)

#But when I am passing it though the CPPFunction it is not working

cppFunction("roll_sd(x, width = 5)")

CodePudding user response：

This is not how cppFunction works. You cannot simply pass it an R call and expect it to magically transpile to C .

In fact, it's more like the other way round. You write the function in C and cppFunction makes that function available in R.

A crude implementation of a rolling standard deviation that matches roll_sd would be something like this:

cppFunction("NumericVector rolling_sd(NumericVector x, int width) {
   NumericVector y = clone(x);
   for(int i = 0; i < (int)x.size(); i  ) {
     if(i < (width - 1)) {
       y[i] = NA_REAL;
     } else {
        double sum = 0.0;
        double total = 0.0;
        for(int j = i - (width - 1); j <= i; j  ) {
          sum  = x[j];
        }
        
        double mean = sum / width;
        
        for(int j = i - (width - 1); j <= i; j  ) {
          total  = pow(x[j] - mean, 2);
        }
        y[i] = sqrt(total / (width - 1));
      }
   }
   return y;
}")

After we run this code, the function rolling_sd is now available to us in R, and gives the same result as roll_sd:

set.seed(1)

x <- rnorm(10)

roll::roll_sd(x, width = 5)
#> [1]        NA        NA        NA        NA 0.9610394 1.0022155 1.0183545
#> [8] 0.8694145 0.6229882 0.6688342

rolling_sd(x, width = 5)
#> [1]        NA        NA        NA        NA 0.9610394 1.0022155 1.0183545
#> [8] 0.8694145 0.6229882 0.6688342

However, our C version is over 10 times faster, as the following benchmark shows.

microbenchmark::microbenchmark(roll_sd(x, width = 5), rolling_sd(x, width = 5))
#> Unit: microseconds
#>                      expr  min   lq   mean median   uq   max neval cld
#>     roll_sd(x, width = 5) 38.7 41.4 44.310   42.2 44.3 154.5   100   b
#>  rolling_sd(x, width = 5)  2.1  2.6  3.346    3.3  3.7  15.0   100  a