Algoritma median bergulir di C

Saat ini saya sedang mengerjakan algoritme untuk menerapkan filter median bergulir (analog dengan filter rata-rata bergulir) di C. Dari penelusuran literatur saya, tampaknya ada dua cara yang cukup efisien untuk melakukannya. Yang pertama adalah mengurutkan jendela nilai awal, kemudian melakukan pencarian biner untuk memasukkan nilai baru dan menghapus nilai yang sudah ada di setiap iterasi.

Yang kedua (dari Hardle dan Steiger, 1995, JRSS-C, Algoritma 296) membangun struktur heap berujung ganda, dengan maxheap di satu ujung, minheap di sisi lain, dan median di tengah. Ini menghasilkan algoritma waktu-linier, bukan yang O (n log n).

Inilah masalah saya: menerapkan yang pertama dapat dilakukan, tetapi saya perlu menjalankannya pada jutaan rangkaian waktu, jadi efisiensi sangat penting. Yang terakhir ini terbukti sangat sulit untuk diterapkan. Saya menemukan kode di file Trunmed.c dari kode untuk paket statistik R, tetapi ini agak tidak dapat diuraikan.

Apakah ada yang tahu tentang implementasi C yang ditulis dengan baik untuk algoritma median waktu linier berguling?

Edit: Tautan ke kode Trunmed.c http://google.com/codesearch/p?hl=id&sa=N&cd=1&ct=rc#mYw3h_Lb_e0/R-2.2.0/src/library/stats/src/Trunmed.c

Saya telah melihat R src/library/stats/src/Trunmed.cbeberapa kali karena saya menginginkan sesuatu yang serupa juga dalam subrutin C ++ class / C mandiri. Perhatikan bahwa ini sebenarnya adalah dua implementasi sekaligus, lihat src/library/stats/man/runmed.Rd(sumber file bantuan) yang menyatakan

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

Akan menyenangkan melihat ini digunakan kembali dengan cara yang lebih mandiri. Apakah Anda menjadi sukarelawan? Saya dapat membantu dengan beberapa bit R.

Sunting 1 : Selain tautan ke versi Trunmed.c yang lebih lama di atas, berikut adalah salinan SVN saat ini

• Srunmed.c (untuk versi Stuetzle)
• Trunmed.c (untuk versi Turlach)
• runmed.R untuk fungsi R yang memanggil ini

Sunting 2 : Ryan Tibshirani memiliki beberapa kode C dan Fortran pada binning median cepat yang mungkin merupakan titik awal yang cocok untuk pendekatan berjendela.

Saya tidak dapat menemukan implementasi modern dari struktur data c ++ dengan statistik pesanan sehingga akhirnya menerapkan kedua ide di tautan pembuat kode teratas yang disarankan oleh MAK ( Editor Pertandingan : gulir ke bawah ke FloatingMedian).

Dua multiset

Ide pertama mempartisi data menjadi dua struktur data (heaps, multisets, dll) dengan O (ln N) per penyisipan / penghapusan tidak memungkinkan kuantil diubah secara dinamis tanpa biaya yang besar. Yaitu kita dapat memiliki median bergulir, atau 75% bergulir tetapi tidak keduanya pada saat yang bersamaan.

Pohon segmen

Ide kedua menggunakan pohon segmen yaitu O (ln N) untuk penyisipan / penghapusan / kueri tetapi lebih fleksibel. Yang terbaik dari semua "N" adalah ukuran rentang data Anda. Jadi jika median bergulir Anda memiliki jendela satu juta item, tetapi datanya bervariasi dari 1..65536, maka hanya 16 operasi yang diperlukan per pergerakan jendela bergulir 1 juta !!

Kode c ++ mirip dengan yang diposting Denis di atas ("Berikut algoritme sederhana untuk data terkuantisasi")

Pohon Statistik Pesanan GNU

Tepat sebelum menyerah, saya menemukan bahwa stdlibc ++ berisi pohon statistik pesanan !!!

Ini memiliki dua operasi penting:

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

Lihat libstdc ++ manual policy_based_data_structures_test (cari "split and join").

Saya telah membungkus pohon untuk digunakan dalam tajuk praktis untuk kompiler yang mendukung typedefs parsial gaya 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

Saya telah melakukan implementasi C di sini . Beberapa detail lebih lanjut ada dalam pertanyaan ini: Rolling median dalam implementasi C - Turlach .

Penggunaan sampel:

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

Saya menggunakan penaksir median tambahan ini:

median += eta * sgn(sample - median)

yang memiliki bentuk yang sama dengan penaksir rata-rata yang lebih umum:

mean += eta * (sample - mean)

Di sini eta adalah parameter kecepatan pembelajaran kecil (misalnya 0.001), dan sgn()merupakan fungsi signum yang mengembalikan salah satu {-1, 0, 1}. (Gunakan konstanta etaseperti ini jika datanya tidak stasioner dan Anda ingin melacak perubahan dari waktu ke waktu; jika tidak, untuk sumber stasioner gunakan sesuatu seperti eta = 1 / nkonvergen, di mana njumlah sampel yang terlihat sejauh ini.)

Juga, saya memodifikasi penaksir median agar berfungsi untuk kuantil sewenang-wenang. Secara umum, fungsi kuantil memberi tahu Anda nilai yang membagi data menjadi dua pecahan: pdan 1 - p. Berikut ini memperkirakan nilai ini secara bertahap:

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

Nilainya pharus di dalam [0, 1]. Ini pada dasarnya menggeser sgn()keluaran simetris fungsi {-1, 0, 1}untuk condong ke satu sisi, mempartisi sampel data menjadi dua nampan berukuran tidak sama (pecahan pdan 1 - pdata masing-masing kurang dari / lebih besar dari perkiraan kuantil). Perhatikan bahwa untuk p = 0.5, ini mengurangi penduga median.

Berikut algoritme sederhana untuk data terkuantisasi (beberapa bulan kemudian):

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

Rolling median dapat ditemukan dengan mempertahankan dua partisi angka.

Untuk memelihara partisi, gunakan Min Heap dan Max Heap.

Max Heap akan berisi angka yang lebih kecil dari sama dengan median.

Min Heap akan berisi angka yang lebih besar dari sama dengan median.

Balancing Constraint: jika jumlah total elemen genap maka kedua heap harus memiliki elemen yang sama.

jika jumlah elemen ganjil maka Max Heap akan memiliki satu elemen lebih banyak dari Min Heap.

Elemen Median: Jika kedua partisi memiliki jumlah elemen yang sama maka median adalah setengah dari jumlah elemen maks dari partisi pertama dan elemen min dari partisi kedua.

Jika tidak, median akan menjadi elemen maks dari partisi pertama.

Algoritma-
1- Ambil dua Heap (1 Min Heap dan 1 Max Heap)
   Max Heap akan berisi setengah dari jumlah elemen
   Min Heap akan berisi jumlah elemen paruh kedua

2- Bandingkan nomor baru dari aliran dengan bagian atas Max Heap, 
   jika lebih kecil atau sama tambahkan angka tersebut di max heap. 
   Jika tidak, tambahkan nomor di Min Heap.

3- jika Min Heap memiliki lebih banyak elemen daripada Max Heap 
   lalu hapus elemen teratas Min Heap dan tambahkan Max Heap.
   jika maks Heap memiliki lebih dari satu elemen daripada di Min Heap 
   lalu hapus elemen teratas dari Max Heap dan tambahkan Min Heap.

4- Jika kedua tumpukan memiliki jumlah elemen yang sama
   median adalah setengah dari jumlah max elemen dari Max Heap dan min elemen dari Min Heap.
   Jika tidak, median akan menjadi elemen maks dari partisi pertama.

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

Mungkin ada baiknya menunjukkan bahwa ada kasus khusus yang memiliki solusi tepat sederhana: ketika semua nilai dalam aliran adalah bilangan bulat dalam rentang yang ditentukan (relatif) kecil. Misalnya, asumsikan semuanya harus berada di antara 0 dan 1023. Dalam kasus ini, cukup tentukan larik 1024 elemen dan hitungan, dan hapus semua nilai ini. Untuk setiap nilai dalam aliran, tambahkan bin yang sesuai dan hitungannya. Setelah aliran berakhir, temukan bin yang berisi hitungan / 2 nilai tertinggi - mudah dilakukan dengan menambahkan bin berurutan mulai dari 0. Dengan menggunakan metode yang sama, nilai urutan peringkat arbitrer dapat ditemukan. (Ada sedikit kerumitan jika perlu mendeteksi saturasi bin dan "meningkatkan" ukuran nampan penyimpanan ke jenis yang lebih besar selama proses.)

Kasus khusus ini mungkin tampak artifisial, tetapi dalam praktiknya sangat umum. Ini juga dapat diterapkan sebagai perkiraan untuk bilangan real jika mereka berada dalam rentang dan tingkat presisi yang "cukup baik" diketahui. Ini akan berlaku untuk hampir semua rangkaian pengukuran pada sekelompok objek "dunia nyata". Misalnya, tinggi atau berat sekelompok orang. Bukan set yang cukup besar? Ini akan bekerja dengan baik untuk panjang atau berat semua (individu) bakteri di planet ini - dengan asumsi seseorang dapat menyediakan data!

Sepertinya saya salah membaca aslinya - yang sepertinya menginginkan median jendela geser, bukan hanya median aliran yang sangat panjang. Pendekatan ini masih berhasil untuk itu. Muat nilai aliran N pertama untuk jendela awal, lalu untuk nilai aliran N + 1, tambahkan bin yang sesuai sambil mengurangi bin yang sesuai dengan nilai aliran ke-0. Dalam hal ini diperlukan untuk mempertahankan nilai N terakhir untuk memungkinkan penurunan, yang dapat dilakukan secara efisien dengan secara siklis menangani larik berukuran N. Karena posisi median hanya dapat berubah sebesar -2, -1,0,1 , 2 di setiap langkah jendela geser, tidak perlu menjumlahkan semua nampan hingga median di setiap langkah, cukup sesuaikan "penunjuk median" tergantung pada nampan sisi mana yang dimodifikasi. Misalnya, jika nilai baru dan nilai yang dihapus berada di bawah median saat ini, maka nilai tersebut tidak berubah (offset = 0). Metode ini rusak ketika N menjadi terlalu besar untuk disimpan dengan nyaman dalam memori.

Jika Anda memiliki kemampuan untuk mereferensikan nilai sebagai fungsi titik waktu, Anda dapat mengambil sampel nilai dengan penggantian, menerapkan bootstrap untuk menghasilkan nilai median yang di-bootstrap dalam interval keyakinan. Ini memungkinkan Anda menghitung perkiraan median dengan efisiensi yang lebih besar daripada terus-menerus menyortir nilai yang masuk ke dalam struktur data.

Bagi yang membutuhkan running median di Java ... PriorityQueue adalah teman Anda. O (log N) masukkan, O (1) median saat ini, dan O (N) hapus. Jika Anda mengetahui distribusi data Anda, Anda dapat melakukan jauh lebih baik daripada ini.

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

Ini adalah salah satu yang dapat digunakan ketika keluaran yang tepat tidak penting (untuk tujuan tampilan, dll.) Anda membutuhkan totalcount dan lastmedian, ditambah nilai baru.

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

Menghasilkan hasil yang cukup tepat untuk hal-hal seperti page_display_time.

Aturan: aliran input harus lancar pada urutan waktu tampilan halaman, jumlah besar (> 30 dll), dan memiliki median bukan nol.

Contoh: waktu buka halaman, 800 item, 10ms ... 3000ms, rata-rata 90ms, median nyata: 11ms

Setelah 30 input, kesalahan median umumnya <= 20% (9ms..12ms), dan semakin berkurang. Setelah 800 masukan, kesalahannya adalah + -2%.

Pemikir lain dengan solusi serupa ada di sini: Median Filter Implementasi super efisien

Berikut adalah implementasi 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
    }
}

Jika Anda hanya membutuhkan rata-rata yang diperhalus, cara cepat / mudah adalah mengalikan nilai terakhir dengan x dan nilai rata-rata dengan (1-x) lalu menjumlahkannya. Ini kemudian menjadi rata-rata baru.

edit: Bukan apa yang diminta pengguna dan tidak valid secara statistik tetapi cukup baik untuk banyak penggunaan.
Saya akan meninggalkannya di sini (terlepas dari suara negatifnya) untuk pencarian!