Algoritmo de la mediana de rodadura en C

Actualmente estoy trabajando en un algoritmo para implementar un filtro de mediana móvil (análogo a un filtro de media rodante) en C. De mi búsqueda de la literatura, parece haber dos formas razonablemente eficientes de hacerlo. El primero es ordenar la ventana de valores inicial, 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 lo primero es factible, pero necesito ejecutar esto en millones de series temporales, por lo que la eficiencia es muy importante. Esto último está demostrando ser muy difícil de implementar. Encontré código en el archivo Trunmed.c del código para el paquete de estadísticas de R, pero es bastante indescifrable.

¿Alguien sabe de una implementación de C bien escrita para el algoritmo de mediana lineal de tiempo rodante?

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 src/library/stats/src/Trunmed.c R src/library/stats/src/Trunmed.c algunas veces, ya que quería algo similar en una subrutina C + clase / C independiente. 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 R bits.

Editar 1 : Además del enlace a la versión anterior de Trunmed.c, aquí están las copias SVN actuales de

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

Edición 2 : Ryan Tibshirani tiene un código C y Fortran en el agrupamiento intermedio 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 estadísticas de orden, así que terminé implementando ambas ideas en el enlace de los mejores codificadores sugerido por MAK ( Match Editorial : desplácese hacia abajo a FloatingMedian).

Dos multiestratos

La primera idea divide los datos en dos estructuras de datos (montones, multisecuencias, etc.) con O (ln N) por inserción / supresión no permite que el cuantil se cambie dinámicamente sin un gran costo. Es decir, podemos tener una mediana progresiva, o un 75% de rodadura, pero no ambas 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, 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, solo se requieren 16 operaciones por movimiento de la ventana móvil de 1 millón.

El código de C ++ es similar al que Denis publicó anteriormente (“Aquí hay un algoritmo simple para datos cuantificados”)

Árboles de estadísticas de orden GNU

Justo antes de darme por vencido, encontré 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 libstdc ++ manual policy_based_data_structures_test (busque “dividir y unir”).

He envuelto el árbol para su uso en un encabezado de conveniencia para los comstackdores compatibles con typedefs de estilo 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 

He hecho una implementación de C aquí . Algunos detalles más están en esta pregunta: Rolling Median en la implementación de C – 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); } } 

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

 """ 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 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, 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 

Yo uso este estimador mediano 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 parámetro de velocidad de aprendizaje pequeño (por ejemplo, 0.001 ), y sgn() es la función de signo que devuelve uno de {-1, 0, 1} . (Use una eta constante como esta si los datos no son estacionarios y desea rastrear los cambios a lo largo del tiempo; de lo contrario, para fonts estacionarias use algo como eta = 1 / n para converger, donde n es el número de muestras vistas hasta ahora. )

Además, modifiqué el estimador mediano para hacerlo funcionar para cuantiles arbitrarios. En general, una función de cuantiles te dice el valor que divide los datos en dos fracciones: p y 1 - p . A continuación, se estima este valor de forma incremental:

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

El valor p debe estar dentro de [0, 1] . Esto esencialmente desplaza la salida simétrica {-1, 0, 1} de la función sgn() para inclinarse hacia un lado, dividiendo las muestras de datos en dos intervalos de tamaños desiguales (las fracciones p y 1 - p de los datos son menores que / mayores que la estimación del cuantil, respectivamente). Tenga en cuenta que para p = 0.5 , esto se reduce al estimador mediano.

La mediana del balanceo se puede encontrar manteniendo dos particiones de números.

Para mantener particiones use Min Heap y Max Heap.

Max Heap contendrá números más pequeños que la mediana.

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

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

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 la misma cantidad de elementos, la mediana será la mitad de la sum 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 de elementos
    Min Heap contendrá la segunda mitad de elementos

 2- Comparar el nuevo número de la secuencia con la parte superior de Max Heap, 
    si es más pequeño o igual, sume ese número en el montón máximo. 
    De lo contrario, agregue el número en Min Heap.

 3- if 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 quite el elemento superior de Max Heap y agregue Min Heap.

 4- Si ambos montones tienen la misma cantidad de elementos, entonces
    la mediana será la mitad de la sum 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 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 tiene la capacidad de referenciar valores en función de los puntos en el tiempo, puede muestrear valores con reemplazo, aplicando bootstrapping para generar un valor mediano bootstrapped dentro de intervalos de confianza. Esto puede permitirle calcular una mediana aproximada con mayor eficiencia que constantemente clasificando los valores entrantes en una estructura de datos.

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

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

Quizás vale la pena señalar que hay un caso especial que tiene una solución exacta simple: cuando todos los valores en la secuencia son enteros dentro de un rango relativamente pequeño. Por ejemplo, suponga que todos deben estar entre 0 y 1023. En este caso, solo defina una matriz de 1024 elementos y un recuento, y borre todos estos valores. Para cada valor en la secuencia, incremente el bin correspondiente y el recuento. Una vez que finaliza la secuencia, encuentre la bandeja que contiene el valor más alto de conteo / 2, que se logra fácilmente al agregar contenedores sucesivos comenzando desde 0. Usando el mismo método, se puede encontrar el valor de una orden de clasificación arbitraria. (Hay una pequeña complicación 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 para números reales si se encuentran dentro de un rango y se conoce un nivel de precisión “lo suficientemente bueno”. Esto sería válido para casi cualquier conjunto de mediciones 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 las longitudes o pesos de todas las bacterias (individuales) en el planeta, ¡suponiendo que alguien pudiera suministrar los datos!

Parece que leí mal el original, lo que parece que quiere una ventana deslizante mediana en lugar de la mediana de un flujo muy largo. Este enfoque todavía funciona para eso. Cargue los primeros N valores de flujo para la ventana inicial, luego para el N + 1º valor de flujo, incremente el intervalo correspondiente mientras disminuye el intervalo correspondiente al 0º valor de flujo. En este caso, es necesario retener los últimos valores N para permitir el decremento, lo que se puede hacer de manera eficiente abordando cíclicamente una matriz de tamaño N. Dado que la posición de la mediana solo puede cambiar por -2, -1,0,1 , 2 en cada paso de la ventana deslizante, no es necesario sumr todas las bandejas hasta la mediana en cada paso, simplemente ajuste el “puntero mediano” dependiendo de las bandejas de los laterales que se hayan modificado. Por ejemplo, si tanto el valor nuevo como el que se elimina caen por debajo de la mediana actual, entonces no cambia (desplazamiento = 0). El método se descompone cuando N se vuelve demasiado grande para mantenerse convenientemente en la memoria.

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

 { 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 uniforme en el orden del tiempo de visualización de la página, grande en conteo (> 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 las medianas generalmente es <= 20% (9ms..12ms) y obtiene menos y menos. Después de 800 entradas, el error es + -2%.

Otro pensador con una solución similar está aquí: Median Filter Implementación super 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 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(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 requiere un promedio suavizado, una forma rápida / fácil es multiplicar el último valor por xy el valor promedio por (1-x) y luego agregarlos. Esto luego se convierte en el nuevo promedio.

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