Aan de slagGa gratis aan de slag

Rollende gemiddelden (in C++)

De gevectoriseerde rollende-gemiddeldefunctie, rollmean3(), presteerde goed maar was erg slecht leesbaar. Het is lastig te zien dat deze code een rollend gemiddelde berekent, wat het lastiger maakt om de functie te debuggen en te onderhouden.

De tweede versie, rollmean2(), was langzamer maar makkelijker te lezen. Als we dit naar C++ vertalen, krijgen we hopelijk zowel leesbaarheid als hoge snelheid.

rollmean2() staat gedefinieerd in je werkruimte; je kunt de definitie printen om jezelf eraan te herinneren hoe het werkt. Je gaat rollmean2() nu naar C++ vertalen en toewijzen aan rollmean4().

Deze oefening maakt deel uit van de cursus

R-code optimaliseren met Rcpp

Cursus bekijken

Oefeninstructies

  • Stel res in als een NumericVector met lengte n en waarden gegeven door de get_na()-methode van NumericVector.
  • Bereken total als de eerste window waarden van x.
  • Bereken het gemiddelde op window - 1 als het totaal gedeeld door de breedte van het venster.
  • Werk in de tweede lus het totaal bij door het i - windowde element van x af te trekken en het ide element van x op te tellen.

Praktische interactieve oefening

Probeer deze oefening eens door deze voorbeeldcode in te vullen.

#include 
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector rollmean4(NumericVector x, int window) {
  int n = x.size();
  
  // Set res as a NumericVector of NAs with length n
  NumericVector res(___, ___::___());
  
  // Sum the first window worth of values of x
  double total = 0.0;
  for(int i = 0; i < window; i++) {
    total += ___;
  }
  
  // Treat the first case seperately
  res[window - 1] = ___ / window;
  
  // Iteratively update the total and recalculate the mean 
  for(int i = window; i < n; i++) {
    // Remove the (i - window)th case, and add the ith case
    total += - ___ + ___;
    // Calculate the mean at the ith position
    res[i] = total / window;
  }
  
  return res;  
}

/*** R
   # Compare rollmean2, rollmean3 and rollmean4   
   set.seed(42)
   x <- rnorm(1e4)
   microbenchmark( 
    rollmean2(x, 4), 
    rollmean3(x, 4), 
    rollmean4(x, 4), 
    times = 5
   )   
*/
Code bewerken en uitvoeren