В настоящее время я работаю над алгоритмом реализации скользящего медианного фильтра (аналог скользящего среднего фильтра) на языке 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
Я несколько раз просматривал 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 <- 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 по fast median binning, который может быть подходящей отправной точкой для оконного подхода.
Я не мог'т найти современная реализация на 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 <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
Я'ве сделали к реализации здесь. Еще несколько подробностей в данный вопрос: 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<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)
Вот эта - это небольшой параметра скорости обучения (например, 0.001
), и функция SGN () - это функция Сигнум, который возвращает одно из
{-1, 0, 1}. (Использовать константу
ета, как это, если данные нестационарные и вы хотите отслеживать изменения с течением времени; в противном случае, для стационарных источников, использовать что-то вроде
ета = 1 / N, чтобы сходятся, где
n` - количество отсчетов видел до сих пор.)
Кроме того, я изменил средний оценщик, чтобы заставить его работать для произвольного квантилей. В общем, функция квантиля представляет собой значение, которое делит данные на две фракции: P
и 1 - п
. Следующие оценки постепенно это значение:
quantile += eta * (sgn(sample - quantile) + 2.0 * p - 1.0)
Значение " Р " должно быть в пределах[0, 1]. Это существенно сдвигает
функция SGN()функция'ы симметричный выход
{-1, 0, 1}склоняться к одной стороне, секционирование данных образцов на две неодинаково размер закрома (фракции
Pи
1 - пданных меньше/больше, чем квантильные оценки, соответственно). Обратите внимание, что для п = 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 кучи.
Медианный элемент: если обе секции имеет одинаковое число элементов, то медианой будет половина суммы максимального элемента из первого раздела и 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<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 элементов и количество, и снимите все эти значения. Для каждого значения в потоке приращение соответствующей bin и граф. После завершения потока находим bin, который содержит количество/2 высшая ценность - легко достигается путем добавления последовательные ячейки, начиная с 0. Используя тот же метод, значение произвольного порядка и ранга Могут быть найдены. (Есть небольшие осложнения в случае обнаружения ОГРН насыщенность и "обновление" в размер бункеры большего типа во время выполнения, потребуется.)
Этот особый случай может показаться искусственной, но на практике это очень часто. Он также может быть применен в качестве приближения для действительных чисел, если они находятся в пределах диапазона и усилитель;"достаточно хорошо" в точности известно. Это позволит проводить практически любые измерения на группы и"и quot реальном мире&; объекты. Например, высоты или веса группы людей. Не достаточно большой набор? Он будет работать так же хорошо для длины или веса всех (отдельных) бактерий на планете - если кто-то может предоставить данные!
Похоже, я неправильно понял оригинал - который выглядит, как он хочет, скользящий средний, а не просто средний, а очень длинный поток. Этот подход по-прежнему работает для этого. Нагрузка в N-й поток значений для начального окна, тогда для N+1-го значения потока приращение соответствующей Бен постоянно уменьшающимся ОГРН, соответствующий значению 0-й поток. Надо в таком случае, чтобы сохранить последние N значений, чтобы уменьшить, что может быть сделано эффективно, циклично обращаясь массив размера N. После установки Медианы можно только изменить на -2,-1,0,1,2 на каждом шагу раздвижные окна, это'т необходимо просуммировать все ячейки до медиану на каждом шагу, просто поменяйте на "средний указатель" в зависимости от того, какая сторона(ы) ячейки были изменены. Например, если новое значение и удаляется опуститься ниже нынешнего среднего уровня, то она не'т изменение (смещение = 0). Метод ломается, когда n становится слишком большим, чтобы удобно держать в памяти.
Для тех, кому нужен бег медиана в Java...PriorityQueue-твой друг. O(зарегистрируйте N) вставка O(1) по текущей медианы, и о(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();
}
}
Если у вас есть возможность ссылаться на значения как на функцию моментов времени, вы можете сделать выборку значений с заменой, применяя bootstrapping для генерации бутстраппированного медианного значения в пределах доверительных интервалов. Это может позволить вам вычислить приближенную медиану с большей эффективностью, чем постоянная сортировка входящих значений в структуре данных.
Вот один, который может быть использован для точного выхода не важна (для целей отображения и т. д.) Вам нужно 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
Вот это реализация 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), затем сложить их. Это и будет новым средним значением.
Редактировать: Не то, что просил пользователь, и не так статистически достоверно, но достаточно хорошо для многих случаев.
Я'оставлю это здесь (несмотря на даунвоты) для поиска!