que - Algoritmo de la mediana de rodadura en C

pedir números hasta que se introduzca uno negativo, y calcular la media.

Question

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

Answer 1

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

Answer 2

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://.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

Answer 3

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

Answer 4

He hecho una implementación de C here . 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);
}
}

Answer 5

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

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

Answer 6

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 compiladores compatibles con typedefs 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

Answer 7

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

Answer 8

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

Answer 9

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!

Answer 10

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.

Answer 11

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 fuentes 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 cuantil (http://en.wikipedia.org/wiki/Quantile_function) le indica 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 contenedores de tamaño desigual (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.

Answer 12

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.