Algorithme médian roulant en C

Je travaille actuellement sur un algorithme pour implémenter un filtre de médiane roulante (analogue à un filtre à moyenne mobile) en C. D’après ma recherche dans la littérature, il semble y avoir deux manières raisonnablement efficaces de le faire. La première consiste à sortinger la fenêtre de valeurs initiale, puis à effectuer une recherche binary pour insérer la nouvelle valeur et supprimer celle existante à chaque itération.

Le second (tiré de Hardle et Steiger, 1995, JRSS-C, algorithme 296) construit une structure de tas à deux extrémités, avec un maxheap à une extrémité, un minheap à l’autre et la médiane au milieu. Cela donne un algorithme de temps linéaire au lieu de celui qui est O (n log n).

Voici mon problème: l’implémentation du premier est faisable, mais je dois l’exécuter sur des millions de séries chronologiques, donc l’efficacité compte beaucoup. Ce dernier s’avère très difficile à mettre en œuvre. J’ai trouvé du code dans le fichier Trunmed.c du code pour le package de statistiques de R, mais il est plutôt indéchiffrable.

Est-ce que quelqu’un connaît une implémentation en C bien écrite pour l’algorithme de la médiane linéaire de roulement de temps?

Edit: Lien vers le code Trunmed.c http://google.com/codesearch/p?hl=fr&sa=N&cd=1&ct=rc#mYw3h_Lb_e0/R-2.2.0/src/library/stats/src/Trunmed.c

J’ai examiné le src/library/stats/src/Trunmed.c plusieurs fois, car je voulais aussi quelque chose de similaire dans une sous-routine C ++ / classe C ++ autonome. Notez qu’il s’agit en fait de deux implémentations en un, voir src/library/stats/man/runmed.Rd (la source du fichier d’aide) qui indique

 \details{ Apart from the end values, the result \code{y = runmed(x, k)} simply has \code{y[j] = median(x[(j-k2):(j+k2)])} (k = 2*k2+1), computed very efficiently. The two algorithms are internally entirely different: \describe{ \item{"Turlach"}{is the Härdle-Steiger algorithm (see Ref.) as implemented by Berwin Turlach. A tree algorithm is used, ensuring performance \eqn{O(n \log k)}{O(n * log(k))} where \code{n < - length(x)} which is asymptotically optimal.} \item{"Stuetzle"}{is the (older) Stuetzle-Friedman implementation which makes use of median \emph{updating} when one observation enters and one leaves the smoothing window. While this performs as \eqn{O(n \times k)}{O(n * k)} which is slower asymptotically, it is considerably faster for small \eqn{k} or \eqn{n}.} } } 

Ce serait bien de voir cela réutilisé de manière plus autonome. Faites-vous du bénévolat? Je peux aider avec certains des bits R.

Edit 1 : Outre le lien vers l'ancienne version de Trunmed.c ci-dessus, voici les copies SVN actuelles de

  • Srunmed.c (pour la version Stuetzle)
  • Trunmed.c (pour la version Turlach)
  • runmed.R pour la fonction R appelant ces

Edit 2 : Ryan Tibshirani a un code C et Fortran sur le binning rapide médian qui peut être un sharepoint départ approprié pour une approche fenêtrée.

Je n’ai pas pu trouver d’implémentation moderne d’une structure de données c ++ avec des statistiques d’ordre, et j’ai donc fini par implémenter les deux idées dans le lien des codeurs de haut niveau suggéré par MAK ( éditorial Match : FloatingMedian).

Deux multisets

La première idée partitionne les données en deux structures de données (tas, multisets, etc.) avec O (lnN) par insertion / suppression ne permet pas de modifier dynamicment le quantile sans coût élevé. C’est-à-dire que nous pouvons avoir une médiane de roulement, ou 75% de roulement, mais pas les deux à la fois.

Arbre de segment

La deuxième idée utilise une arborescence de segment qui est O (lnN) pour les insertions / suppressions / requêtes mais qui est plus flexible. Mieux encore, le “N” correspond à la taille de votre plage de données. Donc, si votre médiane continue a une fenêtre d’un million d’éléments, mais que vos données varient de 1.65536, seules 16 opérations sont requirejses par mouvement de la fenêtre mobile de 1 million!

Le code c ++ est similaire à ce que Denis a posté ci-dessus (“Voici un algorithme simple pour les données quantifiées”)

Arbres de statistiques d’ordre GNU

Juste avant d’abandonner, j’ai trouvé que stdlibc ++ contenait des arbres statistiques de commande !!!

Celles-ci ont deux opérations critiques:

 iter = tree.find_by_order(value) order = tree.order_of_key(value) 

Voir le manuel libstdc ++ policy_based_data_structures_test (recherchez “split and join”).

J’ai encapsulé l’arborescence pour l’utiliser dans un en-tête de commodité pour les compilateurs prenant en charge les typedefs partiels de style c ++ 0x / c ++ 11:

 #if !defined(GNU_ORDER_STATISTIC_SET_H) #define GNU_ORDER_STATISTIC_SET_H #include  #include  // A red-black tree table storing ints and their order // statistics. Note that since the tree uses // tree_order_statistics_node_update as its update policy, then it // includes its methods by_order and order_of_key. template  using t_order_statistic_set = __gnu_pbds::tree< T, __gnu_pbds::null_type, std::less, __gnu_pbds::rb_tree_tag, // This policy updates nodes' metadata for order statistics. __gnu_pbds::tree_order_statistics_node_update>; #endif //GNU_ORDER_STATISTIC_SET_H 

J’ai effectué une implémentation en C ici . Quelques détails supplémentaires se trouvent dans cette question: Médiane mobile dans la mise en œuvre de C-Turlach .

Exemple d’utilisation:

 int main(int argc, char* argv[]) { int i,v; Mediator* m = MediatorNew(15); for (i=0;i<30;i++) { v = rand()&127; printf("Inserting %3d \n",v); MediatorInsert(m,v); v=MediatorMedian(m); printf("Median = %3d.\n\n",v); ShowTree(m); } } 

Voici un algorithme simple pour les données quantifiées (mois plus tard):

 """ median1.py: moving median 1d for quantized, eg 8-bit data Method: cache the median, so that wider windows are faster. The code is simple -- no heaps, no trees. Keywords: median filter, moving median, running median, numpy, scipy See Perreault + Hebert, Median Filtering in Constant Time, 2007, http://nomis80.org/ctmf.html: nice 6-page paper and C code, mainly for 2d images Example: y = medians( x, window=window, nlevel=nlevel ) uses: med = Median1( nlevel, window, counts=np.bincount( x[0:window] )) med.addsub( +, - ) -- see the picture in Perreault m = med.median() -- using cached m, summ How it works: picture nlevel=8, window=3 -- 3 1s in an array of 8 counters: counts: . 1 . . 1 . 1 . sums: 0 1 1 1 2 2 3 3 ^ sums[3] < 2 <= sums[4] <=> median 4 addsub( 0, 1 ) m, summ stay the same addsub( 5, 1 ) slide right addsub( 5, 6 ) slide left Updating `counts` in an `addsub` is sortingvial, updating `sums` is not. But we can cache the previous median `m` and the sum to m `summ`. The less often the median changes, the faster; so fewer levels or *wider* windows are faster. (Like any cache, run time varies a lot, depending on the input.) See also: scipy.signal.medfilt -- runtime roughly ~ window size http://stackoverflow.com/questions/1309263/rolling-median-algorithm-in-c """ from __future__ import division import numpy as np # bincount, pad0 __date__ = "2009-10-27 oct" __author_email__ = "denis-bz-py at t-online dot de" #............................................................................... class Median1: """ moving median 1d for quantized, eg 8-bit data """ def __init__( s, nlevel, window, counts ): s.nlevel = nlevel # >= len(counts) s.window = window # == sum(counts) s.half = (window // 2) + 1 # odd or even s.setcounts( counts ) def median( s ): """ step up or down until sum cnt to m-1 < half <= sum to m """ if s.summ - s.cnt[sm] < s.half <= s.summ: return sm j, sumj = sm, s.summ if sumj <= s.half: while j < s.nlevel - 1: j += 1 sumj += s.cnt[j] # print "j sumj:", j, sumj if sumj - s.cnt[j] < s.half <= sumj: break else: while j > 0: sumj -= s.cnt[j] j -= 1 # print "j sumj:", j, sumj if sumj - s.cnt[j] < s.half <= sumj: break sm, s.summ = j, sumj return sm def addsub( s, add, sub ): s.cnt[add] += 1 s.cnt[sub] -= 1 assert s.cnt[sub] >= 0, (add, sub) if add < = sm: s.summ += 1 if sub <= sm: s.summ -= 1 def setcounts( s, counts ): assert len(counts) <= s.nlevel, (len(counts), s.nlevel) if len(counts) < s.nlevel: counts = pad0__( counts, s.nlevel ) # numpy array / list sumcounts = sum(counts) assert sumcounts == s.window, (sumcounts, s.window) s.cnt = counts s.slowmedian() def slowmedian( s ): j, sumj = -1, 0 while sumj < s.half: j += 1 sumj += s.cnt[j] sm, s.summ = j, sumj def __str__( s ): return ("median %d: " % sm) + \ "".join([ (" ." if c == 0 else "%2d" % c) for c in s.cnt ]) #............................................................................... def medianfilter( x, window, nlevel=256 ): """ moving medians, y[j] = median( x[j:j+window] ) -> a shorter list, len(y) = len(x) - window + 1 """ assert len(x) >= window, (len(x), window) # np.clip( x, 0, nlevel-1, out=x ) # cf http://scipy.org/Cookbook/Rebinning cnt = np.bincount( x[0:window] ) med = Median1( nlevel=nlevel, window=window, counts=cnt ) y = (len(x) - window + 1) * [0] y[0] = med.median() for j in xrange( len(x) - window ): med.addsub( x[j+window], x[j] ) y[j+1] = med.median() return y # list # return np.array( y ) def pad0__( x, tolen ): """ pad x with 0 s, numpy array or list """ n = tolen - len(x) if n > 0: try: x = np.r_[ x, np.zeros( n, dtype=x[0].dtype )] except NameError: x += n * [0] return x #............................................................................... if __name__ == "__main__": Len = 10000 window = 3 nlevel = 256 period = 100 np.set_printoptions( 2, threshold=100, edgeitems=10 ) # print medians( np.arange(3), 3 ) sinwave = (np.sin( 2 * np.pi * np.arange(Len) / period ) + 1) * (nlevel-1) / 2 x = np.asarray( sinwave, int ) print "x:", x for window in ( 3, 31, 63, 127, 255 ): if window > Len: continue print "medianfilter: Len=%d window=%d nlevel=%d:" % (Len, window, nlevel) y = medianfilter( x, window=window, nlevel=nlevel ) print np.array( y ) # end median1.py 

J’utilise cet estimateur incrémental médian:

 median += eta * sgn(sample - median) 

qui a la même forme que l’estimateur moyen le plus commun:

 mean += eta * (sample - mean) 

Ici, eta est un petit paramètre de vitesse d’apprentissage (par exemple, 0.001 ), et sgn() est la fonction signum qui renvoie une des {-1, 0, 1} de {-1, 0, 1} . (Utilisez une constante comme celle-ci si les données ne sont pas stationnaires et que vous souhaitez suivre les changements au fil du temps; sinon, pour les sources fixes, utilisez quelque chose comme eta = 1 / n pour converger, n étant le nombre d’échantillons observés jusqu’à présent. )

De plus, j’ai modifié l’estimateur médian pour le faire fonctionner pour des quantiles arbitraires. En général, une fonction quantile vous indique la valeur qui divise les données en deux fractions: p et 1 - p . Ce qui suit estime cette valeur de manière incrémentielle:

 quantile += eta * (sgn(sample - quantile) + 2.0 * p - 1.0) 

La valeur p doit être comprise entre [0, 1] . Cela déplace essentiellement la sortie symésortingque de la fonction sgn() {-1, 0, 1} pour se pencher d’un côté, en partitionnant les échantillons de données en deux segments de taille inégale (les fractions p et 1 - p des données sont inférieures / supérieures à l’estimation quantile, respectivement). Notez que pour p = 0.5 , cela se réduit à l’estimateur médian.

La médiane de roulement peut être trouvée en maintenant deux partitions de nombres.

Pour maintenir les partitions, utilisez Min Heap et Max Heap.

Max Heap contiendra des nombres inférieurs à la médiane.

Min Heap contiendra des nombres supérieurs à la médiane.

Contrainte d’équilibrage: si le nombre total d’éléments est égal, les deux segments doivent avoir les mêmes éléments.

Si le nombre total d’éléments est impair, Max Heap aura un élément de plus que Min Heap.

Élément médian: Si les deux partitions ont un nombre égal d’éléments, la médiane sera la moitié de la sum de l’élément max de la première partition et de l’élément min de la seconde partition.

Sinon, la médiane sera l’élément max de la première partition.

 Algorithme-
 1- Prenez deux tas (1 Min Heap et 1 Max Heap)
    Max Heap contiendra le premier demi-nombre d'éléments
    Min Heap contiendra la deuxième moitié du nombre d'éléments

 2- Comparer le nouveau numéro du stream avec le haut de Max Heap, 
    si c'est plus petit ou égal, ajoutez ce nombre dans le tas maximum. 
    Sinon, ajoutez le numéro dans Min Heap.

 3- si min Heap a plus d'éléments que Max Heap 
    puis supprimez l'élément supérieur de Min Heap et ajoutez Max Heap.
    si max Heap a plus d'un élément que dans Min Heap 
    puis supprimez l'élément supérieur de Max Heap et ajoutez Min Heap.

 4- Si les deux tas ont le même nombre d'éléments alors
    la médiane sera la moitié de la sum de l'élément max de Max Heap et de l'élément min de Min Heap.
    Sinon, la médiane sera l'élément max de la première partition.
 public class Solution { public static void main(Ssortingng[] args) { Scanner in = new Scanner(System.in); RunningMedianHeaps s = new RunningMedianHeaps(); int n = in.nextInt(); for(int a_i=0; a_i < n; a_i++){ printMedian(s,in.nextInt()); } in.close(); } public static void printMedian(RunningMedianHeaps s, int nextNum){ s.addNumberInHeap(nextNum); System.out.printf("%.1f\n",s.getMedian()); } } class RunningMedianHeaps{ PriorityQueue minHeap = new PriorityQueue(); PriorityQueue maxHeap = new PriorityQueue(Comparator.reverseOrder()); public double getMedian() { int size = minHeap.size() + maxHeap.size(); if(size % 2 == 0) return (maxHeap.peek()+minHeap.peek())/2.0; return maxHeap.peek()*1.0; } private void balanceHeaps() { if(maxHeap.size() < minHeap.size()) { maxHeap.add(minHeap.poll()); } else if(maxHeap.size() > 1+minHeap.size()) { minHeap.add(maxHeap.poll()); } } public void addNumberInHeap(int num) { if(maxHeap.size()==0 || num < = maxHeap.peek()) { maxHeap.add(num); } else { minHeap.add(num); } balanceHeaps(); } } 

Si vous avez la possibilité de référencer des valeurs en fonction de points dans le temps, vous pouvez échantillonner des valeurs avec remplacement, en appliquant un bootstrapping pour générer une valeur médiane amorcée dans les intervalles de confiance. Cela peut vous permettre de calculer une médiane approximative avec plus d’efficacité que de toujours sortinger les valeurs entrantes dans une structure de données.

Pour ceux qui ont besoin d’une médiane en cours d’exécution en Java … PriorityQueue est votre ami. O (log N) insérer, O (1) médiane actuelle et O (N) supprimer. Si vous connaissez la dissortingbution de vos données, vous pouvez faire beaucoup mieux que cela.

 public class RunningMedian { // Two priority queues, one of reversed order. PriorityQueue lower = new PriorityQueue(10, new Comparator() { public int compare(Integer arg0, Integer arg1) { return (arg0 < arg1) ? 1 : arg0 == arg1 ? 0 : -1; } }), higher = new PriorityQueue(); public void insert(Integer n) { if (lower.isEmpty() && higher.isEmpty()) lower.add(n); else { if (n < = lower.peek()) lower.add(n); else higher.add(n); rebalance(); } } void rebalance() { if (lower.size() < higher.size() - 1) lower.add(higher.remove()); else if (higher.size() < lower.size() - 1) higher.add(lower.remove()); } public Integer getMedian() { if (lower.isEmpty() && higher.isEmpty()) return null; else if (lower.size() == higher.size()) return (lower.peek() + higher.peek()) / 2; else return (lower.size() < higher.size()) ? higher.peek() : lower .peek(); } public void remove(Integer n) { if (lower.remove(n) || higher.remove(n)) rebalance(); } } 

Il est peut-être utile de souligner qu’il existe un cas particulier qui a une solution simple et exacte: lorsque toutes les valeurs du stream sont des entiers compris dans une plage (relativement) petite définie. Par exemple, supposons qu’ils doivent tous être compris entre 0 et 1023. Dans ce cas, définissez simplement un tableau de 1024 éléments et un nombre, et effacez toutes ces valeurs. Pour chaque valeur du stream, incrémentez la corbeille correspondante et le compte. Une fois le stream terminé, trouvez la corbeille qui contient la valeur la plus élevée – facilement obtenue en ajoutant des cellules successives à partir de 0. En utilisant la même méthode, la valeur d’un ordre de classement arbitraire peut être trouvée. (Il est nécessaire de détecter la saturation de la corbeille et de mettre à niveau la taille des bacs de stockage au cours d’une parsing si la détection de la saturation des bacs est mineure.)

Ce cas particulier peut sembler artificiel, mais en pratique, il est très courant. On peut également l’appliquer comme une approximation pour des nombres réels s’ils se situent dans une plage et qu’un niveau de précision “suffisant” est connu. Cela permettrait à peu près n’importe quel ensemble de mesures sur un groupe d’objects “du monde réel”. Par exemple, les hauteurs ou les poids d’un groupe de personnes. Pas un ensemble assez grand? Cela fonctionnerait tout aussi bien pour la longueur ou le poids de toutes les bactéries (individuelles) sur la planète – en supposant que quelqu’un puisse fournir les données!

Il semble que j’ai mal interprété l’original – ce qui semble vouloir une fenêtre médiane glissante au lieu de la simple médiane d’un très long stream. Cette approche fonctionne toujours pour cela. Chargez les N premières valeurs de stream pour la fenêtre initiale, puis pour la valeur de stream N + 1, incrémentez la corbeille correspondante tout en décrémentant la corbeille correspondant à la valeur de stream 0. Il est nécessaire dans ce cas de conserver les dernières N valeurs pour permettre la décrémentation, ce qui peut être fait efficacement en adressant cycliquement un tableau de taille N. La position de la médiane ne pouvant changer que de -2, -1,0,1 , 2 à chaque étape de la fenêtre glissante, il n’est pas nécessaire de faire la sum de toutes les cases jusqu’à la médiane à chaque étape, il suffit de régler le “pointeur médian” en fonction du ou des côtés modifiés. Par exemple, si la nouvelle valeur et celle qui sont supprimées sont toutes deux inférieures à la médiane actuelle, elles ne changent pas (offset = 0). La méthode tombe en panne lorsque N devient trop grand pour être stocké facilement en mémoire.

En voici un qui peut être utilisé lorsque la sortie exacte n’est pas importante (à des fins d’affichage, etc.). Vous avez besoin de totalcount et lastmedian, plus la nouvelle valeur.

 { totalcount++; newmedian=lastmedian+(newvalue>lastmedian?1:-1)*(lastmedian==0?newvalue: lastmedian/totalcount*2); } 

Produit des résultats assez exacts pour des choses comme page_display_time.

Règles: le stream d’entrée doit être lisse dans l’ordre d’affichage de la page, avec un grand nombre (> 30, etc.) et une médiane non nulle.

Exemple: temps de chargement de la page, 800 éléments, 10 ms … 3000 ms, moyenne de 90 ms, médiane réelle: 11 ms

Après 30 entrées, l’erreur médiane est généralement < = 20% (9 ms. 12 ms) et devient de moins en moins. Après 800 entrées, l'erreur est de + -2%.

Un autre penseur avec une solution similaire est ici: Median Filter Implémentation super efficace

Voici l’implémentation Java

 package MedianOfIntegerStream; import java.util.Comparator; import java.util.HashSet; import java.util.Iterator; import java.util.Set; import java.util.TreeSet; public class MedianOfIntegerStream { public Set rightMinSet; public Set leftMaxSet; public int numOfElements; public MedianOfIntegerStream() { rightMinSet = new TreeSet(); leftMaxSet = new TreeSet(new DescendingComparator()); numOfElements = 0; } public void addNumberToStream(Integer num) { leftMaxSet.add(num); Iterator iterMax = leftMaxSet.iterator(); Iterator iterMin = rightMinSet.iterator(); int maxEl = iterMax.next(); int minEl = 0; if (iterMin.hasNext()) { minEl = iterMin.next(); } if (numOfElements % 2 == 0) { if (numOfElements == 0) { numOfElements++; return; } else if (maxEl > minEl) { iterMax.remove(); if (minEl != 0) { iterMin.remove(); } leftMaxSet.add(minEl); rightMinSet.add(maxEl); } } else { if (maxEl != 0) { iterMax.remove(); } rightMinSet.add(maxEl); } numOfElements++; } public Double getMedian() { if (numOfElements % 2 != 0) return new Double(leftMaxSet.iterator().next()); else return (leftMaxSet.iterator().next() + rightMinSet.iterator().next()) / 2.0; } private class DescendingComparator implements Comparator { @Override public int compare(Integer o1, Integer o2) { return o2 - o1; } } public static void main(Ssortingng[] args) { MedianOfIntegerStream streamMedian = new MedianOfIntegerStream(); streamMedian.addNumberToStream(1); System.out.println(streamMedian.getMedian()); // should be 1 streamMedian.addNumberToStream(5); streamMedian.addNumberToStream(10); streamMedian.addNumberToStream(12); streamMedian.addNumberToStream(2); System.out.println(streamMedian.getMedian()); // should be 5 streamMedian.addNumberToStream(3); streamMedian.addNumberToStream(8); streamMedian.addNumberToStream(9); System.out.println(streamMedian.getMedian()); // should be 6.5 } } 

Si vous avez simplement besoin d’une moyenne lissée, une méthode simple et rapide consiste à multiplier la dernière valeur par x et la valeur moyenne par (1-x), puis ajoutez-les. Cela devient alors la nouvelle moyenne.

edit: Pas ce que l’utilisateur a demandé et pas aussi statistiquement valide mais assez bon pour beaucoup d’utilisations.
Je vais le laisser ici (malgré les downvotes) pour la recherche!