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

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

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

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

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

Edit: Ссылка на код 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

Комментарии к вопросу (7)

Я несколько раз просматривал R's src/library/stats/src/Trunmed.c, так как хотел получить нечто подобное в отдельном классе C++ / подпрограмме 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 
Комментарии (2)

Я не мог'т найти современная реализация на C++ структуры данных для статистики, так что в итоге реализации идей в топ-кодеров ссылке, предложенной МАК (игре в редакции: прокрутите вниз, чтобы FloatingMedian).

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

Первая мысль разделяет данные на две структуры данных (кучи, мультимножеств и т. д) С О(В Н) на вставку/удаление не позволяет квантиль, чтобы быть динамически изменены без больших затрат. Т. е. мы можем иметь скользящий средний, или переходящий 75%, но не оба одновременно.

Дерево отрезков,

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

Код на C++ похож на то, что Денис выложил выше (на"Здесь'ы простой алгоритм для квантованных данных и")

ГНУ порядке статистические деревьев

Просто прежде чем сдаваться, я обнаружил, что stdlibc++ содержит порядковую статистику деревья!!!

Они имеют две важные операции:

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

См. libstdc++ в ручном policy_based_data_structures_test (поиск и"и на").

Я обернула дерево используемые в заголовке удобство для компиляторов поддерживает 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
Комментарии (2)

Я'ве сделали [к реализации][1][ здесь][2]. Еще несколько подробностей в данный вопрос: https://stackoverflow.com/questions/5527437/rolling-median-in-c-turlach-implementation.

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


int main(int argc, char* argv[])
{
   int i,v;
   Mediator* m = MediatorNew(15);

   for (i=0;i
Комментарии (2)

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

median += eta * sgn(sample - median)

который имеет такую же форму, как более универсальное средство оценки:

mean += eta * (sample - mean)

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

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

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

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

Комментарии (2)

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

""" 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  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
Комментарии (0)

Прокат медиана может быть найден путем поддержания двух разделов цифр.

Для поддержания перегородки использовать кучу мин и Макс кучи.

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

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

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

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

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

В противном случае медиана будет максимальный элемент из первого раздела. в <предварительно> Алгоритм- 1 - Взять две кучи(1 мин 1 Макс кучи и кучи) Максимальный размер кучи будет содержать первую половину количество элементов Мин кучи будет содержать второй половине числа элементов

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

3 - Если мин кучи имеет больше элементов, чем максимальный размер кучи затем удалить верхний элемент из кучи мин и добавить в Максим Кучко. если максимальный размер кучи имеет более чем один элемент, чем в Мин кучи затем удалить верхний элемент максимальный размер кучи и добавьте в Мин кучи.

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 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 
Комментарии (3)

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

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

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

Комментарии (0)

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


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 
Комментарии (4)

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

Комментарии (0)

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

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

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

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

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

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

Еще один мыслитель с подобным решением-это здесь: https://stackoverflow.com/questions/11482529/median-filter-super-efficient-implementation/15150968#15150968

Комментарии (0)

Вот это реализация 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
    }
}
Комментарии (1)

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

Редактировать: Не то, что просил пользователь, и не так статистически достоверно, но достаточно хорошо для многих случаев.
Я'оставлю это здесь (несмотря на даунвоты) для поиска!

Комментарии (3)