Pure R implementation:
psiHampel <- function(r, a, b, c) { psi <- NULL for (i in 1:length(r)) { if (abs(r[i]) <= a) { psi[i] <- r[i] } else if (abs(r[i]) > a & abs(r[i]) <= b) { psi[i] <- a * sign(r[i]) } else if (abs(r[i]) > b & abs(r[i]) <= c) { psi[i] <- (a * (c - abs(r[i]))/(c - b)) * sign(r[i]) } else { psi[i] <- 0 } } return(psi) }
Rcpp/C++ implementation:
#include <Rcpp.h> #include <math.h> using namespace Rcpp; using namespace std; // [[Rcpp::export]] NumericVector psiHampel_rcpp(NumericVector r, double a, double b, double c) { int n = r.size(); NumericVector y = clone(r); for (int i = 0; i < n; i++) { if (abs(r[i]) <= a) { y[i] = r[i]; } else if ( abs(r[i]) > a && abs(r[i]) <= b) { if (r[i] > 0) { y[i] = a; } else { y[i] = -a; } } else if ( abs(r[i]) > b && abs(r[i]) <= c) { if (r[i] > 0) { y[i] = (a * (c - abs(r[i])) / (c - b)); } else { y[i] = (a * (c - abs(r[i])) / (c - b) * -1); } } else { y[i] = 0; } } return y; }
Below is a ggplot2 graph of values a=2, b=4, c=8 where x ranges from -10 to 10:
Both functions have the same result although the C++ implementation is much faster.