Algoritmo de mediana móvil en C

Actualmente estoy trabajando en un algoritmo para implementar un filtro de mediana variable (análogo a un filtro de media variable) en C. A partir de mi búsqueda en la literatura, parece haber dos formas razonablemente eficientes de hacerlo. La primera es ordenar la ventana inicial de valores, luego realizar una búsqueda binaria para insertar el nuevo valor y eliminar el existente en cada iteración.

El segundo (de Hardle y Steiger, 1995, JRSS-C, algoritmo 296) construye una estructura de montón de dos extremos, con un maxheap en un extremo, un minheap en el otro y la mediana en el medio. Esto produce un algoritmo de tiempo lineal en lugar de uno que es O (n log n).

Aquí está mi problema: implementar el primero es factible, pero necesito ejecutarlo en millones de series de tiempo, por lo que la eficiencia importa mucho. Este último está resultando muy difícil de implementar. Encontré código en el archivo Trunmed.c del código del paquete de estadísticas de R, pero es bastante indescifrable.

¿Alguien sabe de una implementación de C bien escrita para el algoritmo de la mediana móvil del tiempo lineal?

Editar: enlace al código de 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

He mirado las R src/library/stats/src/Trunmed.cvarias veces porque también quería algo similar en una subrutina de clase / C independiente de C ++. Tenga en cuenta que en realidad se trata de dos implementaciones en una, consulte src/library/stats/man/runmed.Rd(la fuente del archivo de ayuda) que dice

\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}.}
  }
}

Sería bueno ver esto reutilizado de una manera más independiente. ¿Eres voluntario? Puedo ayudar con algunos de los bits R.

Edición 1 : además del enlace a la versión anterior de Trunmed.c anterior, aquí hay copias actuales de SVN de

• Srunmed.c (para la versión Stuetzle)
• Trunmed.c (para la versión Turlach)
• runmed.R para la función R llamando a estos

Edición 2 : Ryan Tibshirani tiene algo de código C y Fortran en el binning mediano rápido que puede ser un punto de partida adecuado para un enfoque de ventana.

No pude encontrar una implementación moderna de una estructura de datos c ++ con orden-estadística, así que terminé implementando ambas ideas en el enlace de codificadores superiores sugerido por MAK ( Match Editorial : desplácese hacia abajo hasta FloatingMedian).

Dos conjuntos múltiples

La primera idea divide los datos en dos estructuras de datos (montones, conjuntos múltiples, etc.) con O (ln N) por inserción / eliminación no permite que el cuantil se cambie dinámicamente sin un gran costo. Es decir, podemos tener una mediana móvil o un 75% móvil, pero no ambos al mismo tiempo.

Árbol de segmentos

La segunda idea utiliza un árbol de segmentos que es O (ln N) para insertar / eliminar / consultas, pero es más flexible. Lo mejor de todo es que la "N" es el tamaño de su rango de datos. Entonces, si su mediana móvil tiene una ventana de un millón de elementos, pero sus datos varían de 1..65536, ¡entonces solo se requieren 16 operaciones por movimiento de la ventana móvil de 1 millón!

El código c ++ es similar a lo que Denis publicó anteriormente ("Aquí hay un algoritmo simple para datos cuantificados")

Árboles de estadísticas de orden GNU

Justo antes de rendirme, descubrí que stdlibc ++ contiene árboles de estadísticas de orden.

Estos tienen dos operaciones críticas:

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

Consulte el manual de libstdc ++ policy_based_data_structures_test (busque "dividir y unir").

He envuelto el árbol para usarlo en un encabezado de conveniencia para compiladores que admiten definiciones de tipo parciales de estilo 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

He hecho una implementación de C aquí . Algunos detalles más se encuentran en esta pregunta: Mediana móvil en C - Implementación de Turlach .

Uso de muestra:

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

Utilizo este estimador de mediana incremental:

median += eta * sgn(sample - median)

que tiene la misma forma que el estimador de medias más común:

mean += eta * (sample - mean)

Aquí eta es un pequeño parámetro de tasa de aprendizaje (por ejemplo 0.001), y sgn()es la función signum que devuelve uno de {-1, 0, 1}. (Use una constante etacomo esta si los datos no son estacionarios y desea realizar un seguimiento de los cambios a lo largo del tiempo; de lo contrario, para fuentes estacionarias use algo como eta = 1 / nconverger, donden está la cantidad de muestras vistas hasta ahora).

Además, modifiqué el estimador de la mediana para que funcione para cuantiles arbitrarios. En general, una función de cuantiles le dice el valor que divide los datos en dos fracciones: py 1 - p. Lo siguiente estima este valor de forma incremental:

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

El valor pdebe estar dentro [0, 1]. Esto esencialmente desplaza la sgn()salida simétrica de la función {-1, 0, 1}para inclinarse hacia un lado, dividiendo las muestras de datos en dos contenedores de tamaño desigual (las fracciones py 1 - plos datos son menores / mayores que la estimación cuantílica, respectivamente). Tenga en cuenta que para p = 0.5, esto se reduce al estimador mediano.

Aquí hay un algoritmo simple para datos cuantificados (meses después):

""" 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 mediana móvil se puede encontrar manteniendo dos particiones de números.

Para mantener las particiones, use Min Heap y Max Heap.

Max Heap contendrá números menores que iguales a la mediana.

Min Heap contendrá números mayores que iguales a la mediana.

Restricción de equilibrio: si el número total de elementos es par, ambos montones deben tener elementos iguales.

si el número total de elementos es impar, Max Heap tendrá un elemento más que Min Heap.

Elemento mediano: si ambas particiones tienen el mismo número de elementos, la mediana será la mitad de la suma del elemento máximo de la primera partición y el elemento mínimo de la segunda partición.

De lo contrario, la mediana será el elemento máximo de la primera partición.

Algoritmo-
1- Toma dos Heap (1 Min Heap y 1 Max Heap)
   Max Heap contendrá la primera mitad del número de elementos
   Min Heap contendrá la segunda mitad del número de elementos

2- Compare el nuevo número de la transmisión con la parte superior de Max Heap, 
   si es menor o igual, agregue ese número en el montón máximo. 
   De lo contrario, agregue el número en Min Heap.

3- si min Heap tiene más elementos que Max Heap 
   luego elimine el elemento superior de Min Heap y agregue Max Heap.
   si max Heap tiene más de un elemento que en Min Heap 
   luego elimine el elemento superior de Max Heap y agregue Min Heap.

4- Si Ambos montones tienen el mismo número de elementos, entonces
   la mediana será la mitad de la suma del elemento máximo de Max Heap y el elemento mínimo de Min Heap.
   De lo contrario, la mediana será el elemento máximo de la primera partición.

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

Quizás valga la pena señalar que hay un caso especial que tiene una solución exacta simple: cuando todos los valores en la secuencia son números enteros dentro de un rango definido (relativamente) pequeño. Por ejemplo, suponga que todos deben estar entre 0 y 1023. En este caso, simplemente defina una matriz de 1024 elementos y un recuento, y borre todos estos valores. Para cada valor en la secuencia, incremente el contenedor correspondiente y el recuento. Una vez finalizada la transmisión, busque el contenedor que contenga el valor más alto de count / 2, lo que se logra fácilmente agregando contenedores sucesivos comenzando desde 0. Usando el mismo método, se puede encontrar el valor de un orden de rango arbitrario. (Existe una complicación menor si se necesita detectar la saturación del contenedor y "actualizar" el tamaño de los contenedores de almacenamiento a un tipo más grande durante una ejecución).

Este caso especial puede parecer artificial, pero en la práctica es muy común. También se puede aplicar como una aproximación de números reales si se encuentran dentro de un rango y se conoce un nivel de precisión "suficientemente bueno". Esto sería válido para prácticamente cualquier conjunto de medidas en un grupo de objetos del "mundo real". Por ejemplo, las alturas o pesos de un grupo de personas. ¿No es un conjunto lo suficientemente grande? Funcionaría igual de bien para la longitud o el peso de todas las bacterias (individuales) del planeta, ¡asumiendo que alguien pudiera proporcionar los datos!

Parece que leí mal el original, que parece que quiere una mediana de ventana deslizante en lugar de solo la mediana de una secuencia muy larga. Este enfoque todavía funciona para eso. Cargue los primeros N valores de flujo para la ventana inicial, luego, para el valor N + 1 ° flujo, incremente el contenedor correspondiente mientras disminuye el contenedor correspondiente al valor 0 de flujo. En este caso, es necesario retener los últimos valores de N para permitir la disminución, que se puede hacer de manera eficiente direccionando cíclicamente una matriz de tamaño N. Dado que la posición de la mediana solo puede cambiar en -2, -1,0,1 , 2 en cada paso de la ventana deslizante, no es necesario sumar todos los bins hasta la mediana en cada paso, simplemente ajuste el "puntero mediano" dependiendo de qué lado (s) se modificaron los bins. Por ejemplo, si tanto el nuevo valor como el que se está eliminando caen por debajo de la mediana actual, entonces no cambia (compensación = 0). El método falla cuando N se vuelve demasiado grande para guardarlo convenientemente en la memoria.

Si tiene la capacidad de hacer referencia a valores en función de puntos en el tiempo, puede muestrear valores con reemplazo, aplicando bootstrapping para generar un valor medio de bootstrap dentro de intervalos de confianza. Esto puede permitirle calcular una mediana aproximada con mayor eficiencia que ordenar constantemente los valores entrantes en una estructura de datos.

Para aquellos que necesitan una mediana en ejecución en Java ... PriorityQueue es su amigo. O (log N) insertar, O (1) mediana actual y O (N) eliminar. Si conoce la distribución de sus datos, puede hacerlo mucho mejor.

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

Aquí hay uno que se puede usar cuando la salida exacta no es importante (con fines de visualización, etc.). Necesita totalcount y lastmedian, más el newvalue.

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

Produce resultados bastante exactos para cosas como page_display_time.

Reglas: el flujo de entrada debe ser fluido en el orden del tiempo de visualización de la página, contar con un gran recuento (> 30, etc.) y tener una mediana distinta de cero.

Ejemplo: tiempo de carga de la página, 800 elementos, 10 ms ... 3000 ms, promedio de 90 ms, mediana real: 11 ms

Después de 30 entradas, el error de mediana es generalmente <= 20% (9ms..12ms) y es cada vez menor. Después de 800 entradas, el error es + -2%.

Otro pensador con una solución similar está aquí: Median Filter Implementación súper eficiente

Aquí está la implementación de 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 solo necesita un promedio suavizado, una forma rápida / fácil es multiplicar el último valor por x y el valor promedio por (1-x) y luego agregarlos. Este se convierte entonces en el nuevo promedio.

editar: No es lo que pidió el usuario y no es estadísticamente válido, pero lo suficientemente bueno para muchos usos.
¡Lo dejaré aquí (a pesar de los votos negativos) para buscar!