Algoritmo da mediana de rolagem em C

Atualmente, estou trabalhando em um algoritmo para implementar um filtro de média móvel (análogo a um filtro de média móvel) em C. De minha pesquisa na literatura, parece haver duas maneiras razoavelmente eficientes de fazer isso. A primeira é ordenar a janela inicial de valores e, em seguida, realizar uma pesquisa binária para inserir o novo valor e remover o existente a cada iteração.

O segundo (de Hardle e Steiger, 1995, JRSS-C, Algorithm 296) constrói uma estrutura de heap de duas extremidades, com um maxheap em uma extremidade, um minheap na outra e a mediana no meio. Isso produz um algoritmo de tempo linear em vez de um que é O (n log n).

Aqui está o meu problema: implementar o primeiro é possível, mas preciso executá-lo em milhões de séries temporais, portanto, a eficiência é muito importante. Este último está se mostrando muito difícil de implementar. Encontrei o código no arquivo Trunmed.c do código para o pacote de estatísticas do R, mas é indecifrável.

Alguém sabe de uma implementação de C bem escrita para o algoritmo de mediana de rolagem de tempo linear?

Editar: link para o código 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

Eu olhei para R src/library/stats/src/Trunmed.calgumas vezes porque também queria algo semelhante em uma classe C ++ / sub-rotina C autônoma. Observe que, na verdade, são duas implementações em uma, consulte src/library/stats/man/runmed.Rd(a fonte do arquivo de ajuda) que diz

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

Seria bom ver isso reutilizado de uma forma mais autônoma. Você é voluntário? Posso ajudar com alguns dos R bits.

Edição 1 : Além do link para a versão anterior do Trunmed.c acima, aqui estão as cópias atuais do SVN

• Srunmed.c (para a versão Stuetzle)
• Trunmed.c (para a versão Turlach)
• runmed.R para a função R chamando estes

Edição 2 : Ryan Tibshirani tem algum código C e Fortran em binning mediano rápido que pode ser um ponto de partida adequado para uma abordagem em janela.

Não consegui encontrar uma implementação moderna de uma estrutura de dados c ++ com estatística de pedidos, então acabei implementando ambas as ideias no link dos principais codificadores sugerido por MAK ( Match Editorial : role para baixo até FloatingMedian).

Dois multisets

A primeira ideia particiona os dados em duas estruturas de dados (heaps, multisets etc) com O (ln N) por inserção / exclusão não permite que o quantil seja alterado dinamicamente sem um grande custo. Ou seja, podemos ter uma média móvel ou 75% móvel, mas não os dois ao mesmo tempo.

Árvore de segmentos

A segunda ideia usa uma árvore de segmento que é O (ln N) para inserções / exclusões / consultas, mas é mais flexível. O melhor de tudo o "N" é o tamanho do intervalo de dados. Portanto, se sua mediana móvel tem uma janela de um milhão de itens, mas seus dados variam de 1..65536, então apenas 16 operações são necessárias por movimento da janela rolante de 1 milhão !!

O código c ++ é semelhante ao que Denis postou acima ("Aqui está um algoritmo simples para dados quantizados")

Árvores de Estatística de Ordem GNU

Antes de desistir, descobri que stdlibc ++ contém árvores de estatísticas de pedidos !!!

Estes têm duas operações críticas:

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

Veja o manual libstdc ++ policy_based_data_structures_test (procure por "dividir e juntar").

Eu envolvi a árvore para uso em um cabeçalho de conveniência para compiladores que suportam typedefs parciais do 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

Fiz uma implementação C aqui . Mais alguns detalhes estão nesta pergunta: Mediana móvel na implementação C - Turlach .

Uso de amostra:

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

Eu uso este estimador mediano incremental:

median += eta * sgn(sample - median)

que tem a mesma forma que o estimador médio mais comum:

mean += eta * (sample - mean)

Aqui, eta é um pequeno parâmetro de taxa de aprendizagem (por exemplo 0.001), e sgn()é a função signum que retorna um de {-1, 0, 1}. (Use uma constante etacomo esta se os dados não forem estacionários e você quiser rastrear as mudanças ao longo do tempo; caso contrário, para fontes estacionárias, use algo como eta = 1 / nconvergir, onde nestá o número de amostras vistas até agora.)

Além disso, modifiquei o estimador mediano para fazê-lo funcionar para quantis arbitrários. Em geral, uma função de quantil informa o valor que divide os dados em duas frações: pe 1 - p. O seguinte estima esse valor de forma incremental:

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

O valor pdeve estar dentro de [0, 1]. Isso essencialmente muda a sgn()saída simétrica da função {-1, 0, 1}para inclinar para um lado, particionando as amostras de dados em dois compartimentos de tamanhos desiguais (as frações pe 1 - pos dados são menores que / maiores que a estimativa de quantil, respectivamente). Observe que para p = 0.5, isso se reduz ao estimador da mediana.

Aqui está um algoritmo simples para dados quantizados (meses depois):

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

A mediana móvel pode ser encontrada mantendo duas partições de números.

Para manter partições, use Min Heap e Max Heap.

Max Heap conterá números menores que iguais à mediana.

Min Heap conterá números maiores que iguais à mediana.

Restrição de equilíbrio: se o número total de elementos for par, ambos os montes devem ter elementos iguais.

se o número total de elementos for ímpar, o Heap máximo terá um elemento a mais do que o Heap mínimo.

Elemento mediano: se ambas as partições tiverem o mesmo número de elementos, a mediana será a metade da soma do elemento máximo da primeira partição e do elemento mínimo da segunda partição.

Caso contrário, a mediana será o elemento máximo da primeira partição.

Algoritmo-
1- Faça dois Heap (1 Min Heap e 1 Max Heap)
   Max Heap conterá o número de elementos da primeira metade
   Min Heap conterá o número de elementos da segunda metade

2- Compare o novo número do stream com o topo do Max Heap, 
   se for menor ou igual, adicione esse número no heap máximo. 
   Caso contrário, adicione o número em Min Heap.

3- se min Heap tiver mais elementos que Max Heap 
   em seguida, remova o elemento superior de Min Heap e adicione Max Heap.
   se max Heap tiver mais de um elemento do que Min Heap 
   em seguida, remova o elemento superior de Max Heap e adicione Min Heap.

4- Se ambos os heaps tiverem igual número de elementos, então
   a mediana será a metade da soma do elemento máximo do Heap máximo e do elemento mínimo do Heap mínimo.
   Caso contrário, a mediana será o elemento máximo da primeira partição.

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

Talvez valha a pena apontar que existe um caso especial que tem uma solução exata simples: quando todos os valores no fluxo são inteiros dentro de um intervalo definido (relativamente) pequeno. Por exemplo, suponha que todos devem estar entre 0 e 1023. Nesse caso, apenas defina uma matriz de 1024 elementos e uma contagem e apague todos esses valores. Para cada valor no incremento de fluxo, o compartimento correspondente e a contagem. Depois que o fluxo termina, encontre o compartimento que contém o valor mais alto de contagem / 2 - facilmente realizado adicionando recipientes sucessivos a partir de 0. Usando o mesmo método, o valor de uma ordem de classificação arbitrária pode ser encontrado. (Há uma pequena complicação se for necessário detectar a saturação do compartimento e "atualizar" o tamanho dos compartimentos de armazenamento para um tipo maior durante uma execução.)

Este caso especial pode parecer artificial, mas na prática é muito comum. Também pode ser aplicado como uma aproximação para números reais se eles estiverem dentro de um intervalo e um nível "bom o suficiente" de precisão for conhecido. Isso valeria para praticamente qualquer conjunto de medições em um grupo de objetos do "mundo real". Por exemplo, a altura ou o peso de um grupo de pessoas. Não é um conjunto grande o suficiente? Funcionaria igualmente bem para os comprimentos ou pesos de todas as bactérias (individuais) do planeta - supondo que alguém pudesse fornecer os dados!

Parece que eu interpretei mal o original - que parece que ele quer uma mediana de janela deslizante em vez de apenas a mediana de um riacho muito longo. Essa abordagem ainda funciona para isso. Carregue os primeiros N valores de fluxo para a janela inicial e, em seguida, para o N + 1º valor de fluxo, incremente o compartimento correspondente enquanto diminui o compartimento correspondente ao 0º valor de fluxo. É necessário, neste caso, reter os últimos N valores para permitir o decréscimo, o que pode ser feito de forma eficiente endereçando ciclicamente uma matriz de tamanho N. Uma vez que a posição da mediana só pode mudar em -2, -1,0,1 , 2 em cada degrau da janela deslizante, não é necessário somar todos os escaninhos até a mediana de cada degrau, basta ajustar o "ponteiro mediano" dependendo de quais escaninhos laterais foram modificados. Por exemplo, se o novo valor e o que está sendo removido ficarem abaixo da mediana atual, ele não mudará (deslocamento = 0). O método falha quando N se torna muito grande para ser guardado convenientemente na memória.

Se você tiver a capacidade de referenciar valores como uma função de pontos no tempo, poderá amostrar valores com substituição, aplicando bootstrapping para gerar um valor mediano bootstrapped dentro de intervalos de confiança. Isso pode permitir que você calcule uma mediana aproximada com maior eficiência do que classificar constantemente os valores recebidos em uma estrutura de dados.

Para quem precisa de um mediano rodando em Java ... PriorityQueue é seu amigo. Inserção de O (log N), mediana de corrente de O (1) e remoção de O (N). Se você conhece a distribuição de seus dados, pode fazer muito melhor do que isso.

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

Aqui está um que pode ser usado quando a saída exata não é importante (para fins de exibição, etc.). Você precisa de totalcount e lastmedian, mais o novo valor.

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

Produz resultados bastante exatos para coisas como page_display_time.

Regras: o fluxo de entrada precisa ser regular na ordem do tempo de exibição da página, grande em contagem (> 30 etc) e ter uma mediana diferente de zero.

Exemplo: tempo de carregamento da página, 800 itens, 10ms ... 3000ms, média 90ms, mediana real: 11ms

Após 30 entradas, o erro médio é geralmente <= 20% (9ms..12ms) e fica cada vez menor. Após 800 entradas, o erro é + -2%.

Outro pensador com uma solução semelhante está aqui: Median Filter Implementação supereficiente

Aqui está a implementação 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
    }
}

Se você precisar apenas de uma média suavizada, uma maneira rápida / fácil é multiplicar o último valor por xe o valor médio por (1-x) e depois adicioná-los. Isso então se torna a nova média.

editar: Não é o que o usuário pediu e não é estatisticamente válido, mas bom o suficiente para muitos usos.
Vou deixar aqui (apesar dos votos negativos) para pesquisa!