Filtrage des valeurs aberrantes - comment rendre plus rapide la fonction Hampel basée sur la médiane?

Question

Je dois utiliser un filtre Hampel sur mes données, en supprimant les valeurs aberrantes. Je n'ai pas réussi à trouver un existant en Python; que dans Matlab et R.

[description de la fonction Matlab] [1]

[Statistiques discussion Échange de fonction Matlab Hampel] [2]

[R vignette paquet pracma; contient la fonction hampel] [3]

J'ai écrit la fonction suivante, en la modélisant de la fonction dans le paquet R pracma; Cependant, il est beaucoup plus lent que la version Matlab. Ce n'est pas idéal. apprécierait la contribution sur la façon de l'accélérer.

La fonction est indiquée ci-dessous-

def hampel(x,k, t0=3): 
    '''adapted from hampel function in R package pracma 
    x= 1-d numpy array of numbers to be filtered 
    k= number of items in window/2 (# forward and backward wanted to capture in median filter) 
    t0= number of standard deviations to use; 3 is default 
    ''' 
    n = len(x) 
    y = x #y is the corrected series 
    L = 1.4826 
    for i in range((k + 1),(n - k)): 
     if np.isnan(x[(i - k):(i + k+1)]).all(): 
      continue 
     x0 = np.nanmedian(x[(i - k):(i + k+1)]) 
     S0 = L * np.nanmedian(np.abs(x[(i - k):(i + k+1)] - x0)) 
     if (np.abs(x[i] - x0) > t0 * S0): 
      y[i] = x0 
    return(y) 

La mise en œuvre de R dans le paquet « de pracma », que je suis en utilisant comme modèle:

function (x, k, t0 = 3) 
{ 
    n <- length(x) 
    y <- x 
    ind <- c() 
    L <- 1.4826 
    for (i in (k + 1):(n - k)) { 
     x0 <- median(x[(i - k):(i + k)]) 
     S0 <- L * median(abs(x[(i - k):(i + k)] - x0)) 
     if (abs(x[i] - x0) > t0 * S0) { 
      y[i] <- x0 
      ind <- c(ind, i) 
     } 
    } 
    list(y = y, ind = ind) 
} 

Toute aide à rendre plus efficace la fonction, ou un pointeur sur une implémentation existante dans un module Python existant serait très apprécié. Exemple de données ci-dessous %% magie cellulaire timeit en Jupyter indique qu'il faut actuellement 15 secondes pour exécuter:

vals=np.random.randn(250000) 
vals[3000]=100 
vals[200]=-9000 
vals[-300]=8922273 
%%timeit 
hampel(vals, k=6) 

[1]: https://www.mathworks.com/help/signal/ref/hampel.html [2]: https://dsp.stackexchange.com/questions/26552/what-is-a-hampel-filter-and-how-does-it-work [3]: https://cran.r-project.org/web/packages/pracma/pracma.pdf

Answer 1

une solution Pandas est de plusieurs ordres de grandeur plus rapide:

def hampel(vals_orig, k=7, t0=3): 
    ''' 
    vals: pandas series of values from which to remove outliers 
    k: size of window (including the sample; 7 is equal to 3 on either side of value) 
    ''' 
    #Make copy so original not edited 
    vals=vals_orig.copy()  
    #Hampel Filter 
    L= 1.4826 
    rolling_median=vals.rolling(k).median() 
    difference=np.abs(rolling_median-vals) 
    median_abs_deviation=difference.rolling(k).median() 
    threshold= t0 *L * median_abs_deviation 
    outlier_idx=difference>threshold 
    vals[outlier_idx]=np.nan 
    return(vals) 

Temps de synchronisation: 11 ms contre 15 secondes; vaste amélioration.

J'ai trouvé une solution pour un filtre similaire dans this post.