Алгоритм скользящей медианы в C

В настоящее время я работаю над алгоритмом для реализации скользящего медианного фильтра (аналогичного фильтру скользящего среднего) в C. Из моего поиска в литературе, похоже, есть два достаточно эффективных способа сделать это. Первый - отсортировать начальное окно значений, затем выполнить двоичный поиск, чтобы вставить новое значение и удалить существующее на каждой итерации.

Второй (из Hardle and Steiger, 1995, JRSS-C, Algorithm 296) строит двустороннюю структуру кучи с maxheap на одном конце, minheap на другом и медианной средой в середине. Это дает алгоритм линейного времени вместо алгоритма O (n log n).

Вот моя проблема: реализовать первое можно, но мне нужно запустить это на миллионах временных рядов, поэтому эффективность имеет большое значение. Последнее оказывается очень сложно реализовать. Я нашел код в файле Trunmed.c кода для пакета статистики R, но он не поддается расшифровке.

Кто-нибудь знает хорошо написанную реализацию C для алгоритма скользящей медианы с линейным временем?

Изменить: ссылка на код 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

Я смотрел на R src/library/stats/src/Trunmed.cнесколько раз, так как мне тоже хотелось чего-то подобного в автономной подпрограмме C ++ class / C. Обратите внимание, что на самом деле это две реализации в одной, см. src/library/stats/man/runmed.Rd(Источник файла справки), в котором говорится

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

Было бы хорошо, если бы это использовалось повторно в более автономном режиме. Вы работаете волонтером? Я могу помочь с некоторыми R-битами.

Редактировать 1 : Помимо ссылки на старую версию Trunmed.c выше, вот текущие копии SVN

• Srunmed.c (для версии Stuetzle)
• Trunmed.c (для версии Turlach)
• runmed.R для функции R, вызывающей эти

Изменить 2 : у Райана Тибширани есть код C и Fortran для быстрого медианного биннинга, который может быть подходящей отправной точкой для оконного подхода.

Мне не удалось найти современную реализацию структуры данных C ++ со статистикой порядка, поэтому в итоге я реализовал обе идеи в ссылке на топ-кодеры, предложенной MAK ( Match Editor : прокрутите вниз до FloatingMedian).

Два мультимножества

Первая идея разделяет данные на две структуры данных (кучи, мультимножества и т.д.) с O (ln N) на каждую вставку / удаление, что не позволяет динамически изменять квантиль без больших затрат. То есть у нас может быть скользящая медиана или скользящие 75%, но не то и другое одновременно.

Сегментное дерево

Вторая идея использует сегментное дерево O (ln N) для вставки / удаления / запросов, но более гибкое. Лучше всего "N" - это размер вашего диапазона данных. Так что, если ваша скользящая медиана имеет окно из миллиона элементов, но ваши данные варьируются от 1..65536, то для перемещения скользящего окна в 1 миллион требуется всего 16 операций !!

Код на C ++ аналогичен тому, что написал Денис выше («Вот простой алгоритм для квантованных данных»).

Деревья статистики заказов GNU

Незадолго до того, как сдаться, я обнаружил, что stdlibc ++ содержит деревья статистики порядка !!!

У них есть две важные операции:

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

См. Руководство по libstdc ++ policy_based_data_structures_test (поиск по запросу "split and join").

Я обернул дерево для использования в удобном заголовке для компиляторов, поддерживающих частичные определения типов в стиле 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

Я сделал реализацию C здесь . В этом вопросе есть еще несколько деталей: скользящая медиана в реализации C - Turlach .

Пример использования:

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

Я использую эту инкрементальную оценку медианы:

median += eta * sgn(sample - median)

который имеет ту же форму, что и более распространенная оценка среднего:

mean += eta * (sample - mean)

Здесь eta - это небольшой параметр скорости обучения (например, 0.001) и sgn()сигнум-функция, возвращающая одно из значений {-1, 0, 1}. (Используйте такую ​​константу, etaесли данные нестационарны и вы хотите отслеживать изменения с течением времени; в противном случае для стационарных источников используйте что-то вроде eta = 1 / nсходимости, где n- количество отсчетов, просмотренных на данный момент.)

Кроме того, я модифицировал среднюю оценку, чтобы она работала для произвольных квантилей. Как правило, функция квантиля сообщает вам значение, которое делит данные на две части: pи 1 - p. Следующее оценивает это значение постепенно:

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

Значение pдолжно быть в пределах [0, 1]. Это существенно смещает sgn()симметричный вывод функции {-1, 0, 1}в сторону одной стороны, разделяя выборки данных на две ячейки разного размера (доли pи 1 - pданные меньше / больше квантильной оценки, соответственно). Обратите внимание, что для p = 0.5, это сводится к средней оценке.

Вот простой алгоритм для квантованных данных (через несколько месяцев):

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

Скользящую медиану можно найти, поддерживая два разбиения чисел.

Для поддержки разделов используйте Min Heap и Max Heap.

Максимальная куча будет содержать числа, меньшие медианы.

Мин. Куча будет содержать числа, превышающие медианное значение.

Ограничение балансировки: если общее количество элементов четное, то обе кучи должны иметь равные элементы.

если общее количество элементов нечетное, то максимальная куча будет иметь на один элемент больше, чем минимальная куча.

Медианный элемент: если оба раздела имеют одинаковое количество элементов, то медиана будет равна половине суммы максимального элемента из первого раздела и минимального элемента из второго раздела.

В противном случае медиана будет максимальным элементом из первого раздела.

Алгоритм-
1- Возьмите две кучи (1 мин. И 1 макс.)
   Max Heap будет содержать первую половину количества элементов
   Min Heap будет содержать вторую половину количества элементов

2- Сравните новый номер из потока с верхней частью Max Heap, 
   если он меньше или равен, добавьте это число в максимальную кучу. 
   В противном случае добавьте число в Min Heap.

3- если min Heap имеет больше элементов, чем Max Heap 
   затем удалите верхний элемент Min Heap и добавьте Max Heap.
   если max Heap имеет более одного элемента, чем Min Heap 
   затем удалите верхний элемент Max Heap и добавьте Min Heap.

4- Если в обеих кучах одинаковое количество элементов, то
   медиана будет равна половине суммы максимального элемента из максимальной кучи и минимального элемента из минимальной кучи.
   В противном случае медиана будет максимальным элементом из первого раздела.

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

Возможно, стоит указать на особый случай, который имеет простое точное решение: когда все значения в потоке являются целыми числами в пределах (относительно) небольшого определенного диапазона. Например, предположим, что все они должны находиться в диапазоне от 0 до 1023. В этом случае просто определите массив из 1024 элементов и счетчик и очистите все эти значения. Для каждого значения в потоке увеличивают соответствующий интервал и счетчик. После завершения потока найдите ячейку, которая содержит максимальное значение count / 2 - это легко сделать, добавив последовательные ячейки, начиная с 0. Используя тот же метод, можно найти значение произвольного порядка ранжирования. (Есть небольшая сложность, если потребуется обнаружение насыщения бункеров и «обновление» размера бункеров до большего типа во время прогона.)

Этот частный случай может показаться искусственным, но на практике он очень распространен. Его также можно применять в качестве приближения для действительных чисел, если они лежат в пределах диапазона и известен «достаточно хороший» уровень точности. Это справедливо практически для любого набора измерений группы объектов «реального мира». Например, рост или вес группы людей. Недостаточно большой набор? Он будет работать так же хорошо для длины или веса всех (отдельных) бактерий на планете - если кто-то может предоставить данные!

Похоже, я неправильно истолковал оригинал - похоже, что ему нужна медиана скользящего окна, а не просто медиана очень длинного потока. Этот подход по-прежнему работает. Загрузите первые N значений потока для начального окна, затем для значения N + 1-го потока увеличьте соответствующий интервал, уменьшая при этом интервал, соответствующий 0-му значению потока. В этом случае необходимо сохранить последние N значений, чтобы разрешить декремент, что может быть эффективно выполнено путем циклической адресации массива размера N. Поскольку положение медианы может изменяться только на -2, -1,0,1 , 2 на каждом шаге скользящего окна, нет необходимости суммировать все ячейки до медианы на каждом шаге, просто настройте «средний указатель» в зависимости от того, какие стороны (-и) ячейки были изменены. Например, если и новое, и удаляемое значение ниже текущей медианы, оно не меняется (смещение = 0). Метод перестает работать, когда N становится слишком большим, чтобы его можно было удобно хранить в памяти.

Если у вас есть возможность ссылаться на значения в зависимости от моментов времени, вы можете выбрать значения с заменой, применив самозагрузку для генерации начального медианного значения в пределах доверительных интервалов. Это может позволить вам вычислить приблизительную медиану с большей эффективностью, чем постоянная сортировка входящих значений в структуру данных.

Для тех, кому нужна работающая медиана на Java ... PriorityQueue - ваш друг. O (log N) вставить, O (1) текущая медиана и O (N) удалить. Если вы знаете распределение ваших данных, вы можете сделать это намного лучше.

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

Вот тот, который можно использовать, когда точный вывод не важен (для целей отображения и т. Д.). Вам нужны totalcount и lastmedian, а также новое значение.

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

Дает довольно точные результаты для таких вещей, как page_display_time.

Правила: входной поток должен быть плавным в соответствии с порядком времени отображения страницы, большим по счетчику (> 30 и т. Д.) И иметь ненулевую медиану.

Пример: время загрузки страницы, 800 элементов, 10 мс ... 3000 мс, в среднем 90 мс, реальная медиана: 11 мс

После 30 вводов ошибка медианы обычно составляет <= 20% (9 мс..12 мс) и становится все меньше и меньше. После 800 входов погрешность составляет + -2%.

Другой мыслитель с похожим решением находится здесь: Медианный фильтр Суперэффективная реализация

Вот реализация 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
    }
}

Если вам просто требуется сглаженное среднее, быстрый / простой способ - умножить последнее значение на x, а среднее значение на (1-x), а затем сложить их. Затем это становится новым средним значением.

edit: Не то, о чем просил пользователь, и не так статистически достоверно, но достаточно хорошо для множества применений.
Я оставлю его здесь (несмотря на отрицательные голоса) для поиска!