Adding Hampel Psi Function to robustreg

I am currently working on adding the Hampel psi function to robustreg.





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.

No comments:

Post a Comment