Algorithme de la médiane mobile en C

j'ai regardé R's src/library/stats/src/Trunmed.c plusieurs fois car je voulais quelque chose de similaire aussi 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 voir ce produit réutilisé de façon plus autonome. Êtes-vous du bénévolat? Je peux vous aider avec les r bits.

Edition 1 : en plus de la lien vers l'ancienne version de Trunmed.c ci-dessus, voici des copies SVN actuelles de

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

Edit 2 : Ryan Tibshirani a quelques C et Code Fortran sur fast médiane binning , qui peut être un point de départ approprié pour une fenêtre approche.

Je n'ai pas pu trouver une implémentation moderne d'une structure de données c++ avec des statistiques d'ordre donc j'ai fini par implémenter les deux idées dans le lien des codeurs top suggéré par MAK ( Match Editorial : faites défiler Jusqu'à FloatingMedian).

deux multisets

la première idée divise 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 grand coût. I. e. nous pouvons avoir un rouleau médian, ou un rouleau à 75%, mais pas les deux en même temps.

Segment de l'arbre

la deuxième idée utilise une arborescence de segments qui est O(ln N) pour insert/deletions/queries mais est plus flexible. Le meilleur de tous les "N" est la taille de votre plage de données. Donc, si votre rouleau médian a une fenêtre d'un million d'articles, mais vos données varie de 1..65536, alors seulement 16 opérations sont nécessaires par mouvement de la fenêtre de laminage de 1 millions!!

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 L'ordre GNU

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

ils ont deux opérations critiques:

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

Voir libstdc++ manuel policy_based_data_structures_test (recherche de "split et join").

j'ai enveloppé l'arbre pour l'utiliser dans un en-tête de commodité pour les compilateurs supportant C++0x/c++11 style partial typedefs:

#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 un C mise en œuvre ici . Un peu plus de détails dans cette question: Rolling médian en C - Turlach mise en œuvre .

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, 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
    /q/rolling-median-algorithm-in-c-22299/"""

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

j'utilise cet estimateur différentiel de la médiane:

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 eta comme ceci si les données ne sont pas stationnaires et que vous voulez suivre les changements dans le temps; sinon, pour les sources fixes, utilisez quelque chose comme eta = 1 / n pour converger, où n est le nombre d'échantillons observés à ce jour.)

aussi, j'ai modifié l'estimateur médian pour le faire fonctionner pour les quantiles arbitraires. En général, une fonction de quantile vous indique la valeur qui divise les données en deux fractions: p et 1 - p . Les estimations suivantes cette valeur progressivement:

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

la valeur p doit être compris dans [0, 1] . Cela déplace essentiellement la sortie symétrique de la fonction sgn() {-1, 0, 1} vers un côté, en divisant les échantillons de données en deux bacs de taille inégale (les fractions p et 1 - p des données sont inférieures ou supérieures à l'estimation de quantile, respectivement). Notez que pour p = 0.5 , cela se réduit à l'estimateur médian.

la médiane mobile peut être trouvée en maintenant deux partitions de nombres.

pour l'entretien des cloisons utiliser Min Heap et Max Heap.

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

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

Équilibrage De La Contrainte: si le nombre total d'éléments sont même les deux tas devraient avoir les mêmes éléments.

si le nombre total d'éléments est impair alors Max Heap aura un élément de plus que Min Heap.

élément médian: si les deux partitions ont le même nombre d'éléments, alors 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 max de la première partition.

Algorithm-
1- Take two Heap(1 Min Heap and 1 Max Heap)
   Max Heap will contain first half number of elements
   Min Heap will contain second half number of elements

2- Compare new number from stream with top of Max Heap, 
   if it is smaller or equal add that number in max heap. 
   Otherwise add number in Min Heap.

3- if min Heap has more elements than Max Heap 
   then remove top element of Min Heap and add in Max Heap.
   if max Heap has more than one element than in Min Heap 
   then remove top element of Max Heap and add in Min Heap.

4- If Both heaps has equal number of elements then
   median will be half of sum of max element from Max Heap and min element from Min Heap.
   Otherwise median will be max element from the first 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();
    }
}

si vous avez la capacité de faire référence à des valeurs en fonction de points dans le temps, vous pouvez échantillonner des valeurs avec remplacement, en appliquant bootstrapping pour générer une valeur médiane bootstrapped dans les intervalles de confiance. Cela peut vous permettre de calculer une médiane approximative avec une plus grande efficacité que le tri constant des 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 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();
  }
}

il est peut-être intéressant de souligner qu'il existe un cas spécial qui a une solution exacte simple: quand toutes les valeurs dans le flux sont des entiers dans une (relativement) petite gamme définie. Par exemple, supposons qu'ils doivent tous se situent entre 0 et 1023. Dans ce cas, il suffit de définir un tableau de 1024 éléments et un nombre, et effacer toutes ces valeurs. Pour chaque valeur dans le flux, incrémentez la bin correspondante et le nombre. Après la fin du flux trouver la bin qui contient le compte / 2 valeur la plus élevée - facilement obtenue en ajoutant des bacs successifs à partir de 0. En utilisant la même méthode, la valeur d'un ordre de rang arbitraire peut être trouvée. (Il y a une petite complication si la détection de la saturation des bacs et la "mise à niveau" de la taille des bacs de stockage à un plus grand type pendant un passage seront nécessaires.)

Ce cas particulier peut sembler artificiel, mais dans la pratique, il est très commun. Il peut également être appliqué comme une approximation pour les nombres réels s'ils se trouvent dans une gamme et un niveau de précision "suffisant" 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, la hauteur ou le poids d'un groupe de personnes. Pas un assez grand ensemble? Cela fonctionnerait aussi bien pour les longueurs ou les poids de toutes les bactéries (individuelles) sur la planète - en supposant que quelqu'un pourrait fournir les données!

il semble que j'ai mal lu l'original - qui semble comme il veut une médiane de fenêtre coulissante au lieu de la juste la médiane d'un très long flux. Cette approche fonctionne toujours pour cela. Chargez les premières valeurs de n stream pour la fenêtre initiale, puis pour la valeur N+1th stream incrémenter la bin correspondante tout en décrémentant la bin correspondant à la valeur 0th stream. Il est nécessaire dans ce cas de conserver les dernières valeurs de N pour permettre le décrément, ce qui peut être fait efficacement en s'adressant de façon cyclique à un tableau de taille N. puisque la position de la médiane ne peut changer que de -2, -1,0,1,2 sur chaque pas de la fenêtre coulissante, il n'est pas nécessaire de faire la somme de toutes les cases jusqu'à la médiane sur chaque pas, il suffit d'ajuster le "pointeur médian" en fonction de quel côté(s) les cases ont été modifiées. Par exemple, si la nouvelle valeur et celle qui est supprimée tombent sous la médiane courante, alors elle ne change pas (offset = 0). La méthode se décompose Quand N devient trop grand pour tenir commodément dans la mémoire.

En voici un qui peut être utilisé lorsque la sortie exacte n'est pas importante (pour des raisons 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 flux d'entrée doit être lisse sur 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, 10ms...3000ms, moyenne 90ms, médiane réelle: 11ms

après 30 entrées, l'erreur médianes est généralement <=20% (9ms..12ms), et obtient 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 Super efficient implementation

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 est de 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 comme statistiquement valable mais assez bon pour beaucoup d'utilisations.

Je vais le laisser ici (malgré les mauvaises notes) pour la recherche!