Algorithme de médiane glissante en C

Je travaille actuellement sur un algorithme pour implémenter un filtre médian roulant (analogue à un filtre à moyenne mobile) en C. D'après ma recherche dans la littérature, il semble y avoir deux façons raisonnablement efficaces de le faire. La première consiste à trier la fenêtre initiale de valeurs, puis à effectuer une recherche binaire pour insérer la nouvelle valeur et supprimer l'existante à chaque itération.

Le second (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 en temps linéaire au lieu d'un algorithme O (n log n).

Voici mon problème: mettre en œuvre le 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 du package stats de R, mais il est plutôt indéchiffrable.

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

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

J'ai regardé les R à src/library/stats/src/Trunmed.cquelques reprises car je voulais aussi quelque chose de similaire dans un sous-programme autonome de classe / C C ++. Notez qu'il s'agit en fait de deux implémentations en une, voir src/library/stats/man/runmed.Rd(la source du fichier d'aide) qui dit

\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 le voir réutilisé de manière plus autonome. Faites-vous du bénévolat? Je peux vous 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 du code C et Fortran sur le binning médian rapide qui peut être un point de départ approprié pour une approche fenêtrée.

Je n'ai pas pu trouver une 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 meilleurs codeurs suggéré par MAK ( Match Editorial : faites défiler jusqu'à FloatingMedian).

Deux multisets

La première idée partitionne les données en deux structures de données (tas, multisets, etc.) avec O (ln N) par insertion / suppression ne permet pas de modifier dynamiquement le quantile sans un coût élevé. C'est-à-dire que nous pouvons avoir une médiane glissante, ou 75% glissants, mais pas les deux en même temps.

Arborescence des segments

La deuxième idée utilise un arbre de segments qui est O (ln N) pour les insertions / suppressions / requêtes mais qui est plus flexible. Mieux encore, le "N" est la taille de votre plage de données. Donc, si votre médiane glissante a une fenêtre d'un million d'éléments, mais que vos données varient de 1..65536, alors seulement 16 opérations sont nécessaires par mouvement de la fenêtre glissante 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 statistiques de commande GNU

Juste avant d'abandonner, j'ai trouvé que stdlibc ++ contient des arbres de statistiques d'ordre !!!

Ceux-ci ont deux opérations critiques:

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

Voir le manuel de libstdc ++ policy_based_data_structures_test (recherchez "split and join").

J'ai enveloppé l'arborescence pour l'utiliser dans un en-tête pratique 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 <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

// 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 <typename T>
using t_order_statistic_set = __gnu_pbds::tree<
                                  T,
                                  __gnu_pbds::null_type,
                                  std::less<T>,
                                  __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 fait une implémentation C ici . Quelques détails supplémentaires sont dans cette question: Médiane mobile dans l'implémentation 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);
   }
}

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

median += eta * sgn(sample - median)

qui a la même forme que l'estimateur moyen le plus courant:

mean += eta * (sample - mean)

Ici, eta est un petit paramètre de taux d'apprentissage (par exemple 0.001), et sgn()est la fonction signum qui renvoie l'un des {-1, 0, 1}. (Utilisez une constante etacomme 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 stationnaires, utilisez quelque chose comme eta = 1 / npour converger, oùn est le nombre d'échantillons vus jusqu'à présent.)

De plus, j'ai modifié l'estimateur médian pour qu'il fonctionne pour des quantiles arbitraires. En général, une fonction quantile vous indique la valeur qui divise les données en deux fractions: pet 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 pdoit être à l'intérieur [0, 1]. Cela déplace essentiellement la sgn()sortie symétrique de la fonction {-1, 0, 1}vers un côté, en partitionnant les échantillons de données en deux bacs de taille inégale (les fractions pet 1 - pdes 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.

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

""" median1.py: moving median 1d for quantized, e.g. 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 trivial, 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, e.g. 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[s.m] < s.half <= s.summ:
            return s.m
        j, sumj = s.m, 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
        s.m, s.summ = j, sumj
        return s.m

    def addsub( s, add, sub ):
        s.cnt[add] += 1
        s.cnt[sub] -= 1
        assert s.cnt[sub] >= 0, (add, sub)
        if add <= s.m:
            s.summ += 1
        if sub <= s.m:
            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]
        s.m, s.summ = j, sumj

    def __str__( s ):
        return ("median %d: " % s.m) + \
            "".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

La médiane mobile peut être trouvée en conservant 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 pair, les deux tas doivent avoir des éléments égaux.

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 somme de l'élément max de la première partition et de l'élément min de la deuxième partition.

Sinon, la médiane sera l'élément maximum de la première partition.

Algorithme-
1- Prenez deux tas (1 tas min et 1 tas max)
   Max Heap contiendra la première moitié du nombre d'éléments
   Min Heap contiendra la deuxième moitié du nombre d'éléments

2- Comparez le nouveau numéro du flux avec le haut de Max Heap, 
   s'il est plus petit ou égal, ajoutez ce nombre dans le tas max. 
   Sinon, ajoutez un nombre dans Min Heap.

3- si le tas min a plus d'éléments que le tas max 
   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 somme 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 maximum de la première partition.

public class Solution {

    public static void main(String[] 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<Integer> minHeap = new PriorityQueue<Integer>();
    PriorityQueue<Integer> maxHeap = new PriorityQueue<Integer>(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();
    }
}

Il est peut-être intéressant de souligner qu'il existe un cas particulier qui a une solution exacte simple: lorsque toutes les valeurs du flux sont des entiers dans une plage définie (relativement) petite. 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 flux, incrémentez le bac correspondant et le nombre. Une fois le flux terminé, trouvez le bac qui contient la valeur count / 2 la plus élevée - facilement accompli en ajoutant des bacs successifs à partir de 0. En utilisant la même méthode, la valeur d'un ordre de classement arbitraire peut être trouvée. (Il y a une complication mineure si la détection de la saturation du bac et la «mise à niveau» de la taille des bacs de stockage vers un type plus grand pendant une analyse sont nécessaires.)

Ce cas particulier peut sembler artificiel, mais en pratique, il est très courant. Il peut également être appliqué comme une approximation pour les nombres réels s'ils se situent dans une plage et qu'un niveau de précision «assez bon» est connu. Cela vaut pour à peu près n'importe quel ensemble de mesures sur un groupe d'objets «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) de la planète - en supposant que quelqu'un puisse fournir les données!

Il semble que j'ai mal lu l'original - ce qui semble vouloir une médiane de fenêtre glissante au lieu de la médiane d'un très long flux. Cette approche fonctionne toujours pour cela. Chargez les N premières valeurs de flux pour la fenêtre initiale, puis pour la valeur de N + 1ème flux, incrémentez le bac correspondant tout en décrémentant le bac correspondant à la 0ème valeur de flux. Il faut dans ce cas retenir les N dernières valeurs pour permettre la décrémentation, ce qui peut être fait efficacement en adressant cycliquement un tableau de taille N. Puisque la position de la médiane ne peut changer que de -2, -1,0,1 , 2 à chaque étape de la fenêtre glissante, il n'est pas nécessaire de faire la somme de toutes les cases jusqu'à la médiane à chaque étape, il suffit d'ajuster le "pointeur médian" en fonction du (des) côté (s) des cases modifiées. Par exemple, si la nouvelle valeur et celle supprimée tombent en dessous de la médiane actuelle, cela ne change pas (offset = 0). La méthode échoue lorsque N devient trop grand pour être conservé en mémoire.

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 le bootstrap pour générer une valeur médiane bootstrap dans les intervalles de confiance. Cela peut vous permettre de calculer une médiane approximative avec une plus grande efficacité que de trier constamment 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. Insertion O (log N), médiane actuelle O (1) et suppression O (N). Si vous connaissez la distribution de vos données, vous pouvez faire beaucoup mieux que cela.

public class RunningMedian {
  // Two priority queues, one of reversed order.
  PriorityQueue<Integer> lower = new PriorityQueue<Integer>(10,
          new Comparator<Integer>() {
              public int compare(Integer arg0, Integer arg1) {
                  return (arg0 < arg1) ? 1 : arg0 == arg1 ? 0 : -1;
              }
          }), higher = new PriorityQueue<Integer>();

  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();
  }
}

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 newvalue.

{
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 flux d'entrée doit être fluide dans l'ordre du temps d'affichage de la page, grand en nombre (> 30, etc.), et avoir une médiane non nulle.

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

Après 30 entrées, l'erreur médiane est généralement <= 20% (9ms..12ms), 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: Filtre médian Mise en œuvre 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<Integer> rightMinSet;
    public Set<Integer> leftMaxSet;
    public int numOfElements;

    public MedianOfIntegerStream() {
        rightMinSet = new TreeSet<Integer>();
        leftMaxSet = new TreeSet<Integer>(new DescendingComparator());
        numOfElements = 0;
    }

    public void addNumberToStream(Integer num) {
        leftMaxSet.add(num);

        Iterator<Integer> iterMax = leftMaxSet.iterator();
        Iterator<Integer> 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<Integer> {
        @Override
        public int compare(Integer o1, Integer o2) {
            return o2 - o1;
        }
    }

    public static void main(String[] 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 juste besoin d'une moyenne lissée, un moyen rapide / facile consiste à multiplier la dernière valeur par x et la valeur moyenne par (1-x), puis les ajouter. Cela devient alors la nouvelle moyenne.

edit: Pas ce que l'utilisateur a demandé et pas aussi statistiquement valide mais assez bon pour de nombreuses utilisations.
Je vais le laisser ici (malgré les votes négatifs) pour la recherche!