c - standard - r std




Algoritmo mediano de rolagem em C (8)

Atualmente, estou trabalhando em um algoritmo para implementar um filtro mediano de rolagem (análogo a um filtro de rolagem) em C. De minha pesquisa na literatura, parece haver duas maneiras razoavelmente eficientes de fazê-lo. A primeira é classificar a janela inicial de valores e, em seguida, executar uma pesquisa binária para inserir o novo valor e remover o existente em cada iteração.

O segundo (de Hardle e Steiger, 1995, JRSS-C, Algorithm 296) constrói uma estrutura de heap com duas extremidades, com um máximo de um lado, um mínimo sobre o outro 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 é factível, mas eu preciso rodar isso em milhões de séries temporais, então a eficiência é muito importante. Este último está sendo muito difícil de implementar. Eu 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 C bem escrita para o algoritmo mediano de rolling linear?

Editar: vincular ao 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


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

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

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

Min Heap conterá números maiores que igual a mediana.

Restrição de Balanceamento: se o número total de elementos for par, ambos devem ter elementos iguais.

Se o número total de elementos for ímpar, o Max Heap terá mais um elemento que Min Heap.

Elemento Mediano: Se Ambas as partições tiverem o mesmo número de elementos, a mediana será a metade da soma do elemento max da primeira partição e o elemento min da segunda partição.

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

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

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

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, além do novo valor.

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

Produz resultados bastante precisos para coisas como page_display_time.

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

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

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

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


Eu fiz uma implementação de C here . Mais alguns detalhes estão nesta pergunta: Rolling median in C - Turlach implementation .

Uso da 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 aprendizado (por exemplo, 0.001 ), e sgn() é a função signum que retorna um de {-1, 0, 1} . (Use uma constante eta como essa se os dados forem não-estacionários e você quiser acompanhar as mudanças ao longo do tempo; caso contrário, para fontes estacionárias use algo como eta = 1 / n para convergir, onde n é o número de amostras vistas até o momento. )

Além disso, modifiquei o estimador mediano para fazê-lo funcionar para quantis arbitrários. Em geral, uma função quantil informa o valor que divide os dados em duas frações: p e 1 - p . As seguintes estimativas deste valor incrementalmente:

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

O valor p deve estar dentro de [0, 1] . Isso essencialmente desloca a saída simétrica da função sgn() {-1, 0, 1} para se inclinar para um lado, particionando as amostras de dados em duas caixas de tamanho desiguais (frações p e 1 - p dos dados são menores que / maiores que a estimativa do quantil, respectivamente). Note que para p = 0.5 , isso reduz para o estimador mediano.


Para quem precisa de uma mediana de corrida em Java ... PriorityQueue é seu amigo. O (log N) insert, O (1) mediana atual, e O (N) remove. Se você conhece a distribuição dos seus dados, pode fazer muito melhor 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();
  }
}

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


Talvez valha a pena ressaltar que há 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 eles devem estar todos entre 0 e 1023. Nesse caso, apenas defina uma matriz de 1024 elementos e uma contagem e limpe todos esses valores. Para cada valor no fluxo, incremente a caixa correspondente e a contagem. Depois que o fluxo terminar, encontre o bin que contém o maior valor de count / 2 - facilmente obtido adicionando escaninhos 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 detectar a saturação do escaninho e "atualizar" o tamanho dos escaninhos de armazenamento para um tipo maior durante uma execução será necessário.)

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

Parece que eu interpretei mal o original - o que parece que quer uma mediana de janela deslizante em vez de apenas a mediana de um fluxo muito longo. Essa abordagem ainda funciona para isso. Carregue os primeiros N valores de fluxo para a janela inicial e, em seguida, para o valor do fluxo de N + 1ª, incremente a caixa correspondente enquanto diminui a posição correspondente ao valor do fluxo de 0º. É necessário, neste caso, reter os últimos N valores para permitir o decréscimo, o que pode ser feito eficientemente, 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 etapa da janela deslizante, não é necessário somar todos os compartimentos até a mediana de cada etapa, basta ajustar o "ponteiro mediano" dependendo de quais dos lados foram modificados. Por exemplo, se o novo valor e o que está sendo removido ficarem abaixo da mediana atual, ele não mudará (offset = 0). O método é quebrado quando N se torna muito grande para ser armazenado convenientemente na memória.





median