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:
z=seq(-50,50,.01)
all.equal(psiHampel(z,2,4,8), psiHampel_rcpp(z,2,4,8))
[1] TRUE
Both functions have the same result although the C++ implementation is much faster.