我目前正在编写一个用于在 C 中实现滚动中位数滤波器（类似于滚动平均滤波器）的算法。通过查阅文献，似乎有两种比较有效的方法来实现它。第一种方法是对初始值窗口进行排序，然后执行二分搜索以在每次迭代时插入新值并删除现有值。
第二种方法（来自 Hardle 和 Steiger，在 1995 年的 JRSS-C，算法 296 中）构建了一个双端堆结构，其中一个端点是最大堆，另一个端点是最小堆，并且中位数位于中间。这会产生一个线性时间算法，而不是 O(n log n) 的算法。
我的问题在于：虽然第一种方法可行，但我需要在数百万个时间序列上运行此程序，因此效率非常重要。而实现第二种方法却很困难。我在 R 统计软件包的 stats 包中的 Trunmed.c 文件中找到了代码，但它相当难以理解。
有人知道一个良好编写的 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++ 类 / 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 版本:
编辑2: Ryan Tibshirani 在 快速中位数分箱 中有一些C和Fortran代码，可能是窗口方法的合适起点。
我找不到一种带有顺序统计的现代化实现的c++数据结构，因此最终在MAK建议的top coders链接中实现了两个想法（比赛编者按：向下滚动到FloatingMedian）。
第一个想法将数据分成两个数据结构（堆、multiset等），每次插入/删除的复杂度为O(ln N)，但是不允许动态更改分位数而不付出巨大代价。也就是说，我们可以有一个滚动的中位数或滚动的75%，但不能同时有。
第二个想法使用区间树，插入/删除/查询的复杂度为O(ln N)，但更加灵活。最好的是，“N”是您的数据范围的大小。因此，如果您的滚动中位数有100万个项目的窗口，但您的数据变化范围为1..65536，则每移动100万个滚动窗口只需要执行16个操作！
c++代码类似于Denis上面发布的代码（“这是一种用于量化数据的简单算法”）。
就在放弃之前，我发现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风格的部分typedef的编译器使用：
#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实现。更多细节请参见此问题：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)
这里的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
https://dev59.com/73M_5IYBdhLWcg3wmkUK
"""
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
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();
}
}
对于那些需要在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();
}
}
}), higher = new PriorityQueue<Integer>(); 或者 new PriorityQueue<Integer>(10,。我无法运行代码。 - Hengameh这是一个可以在精确输出不重要时使用的函数（用于显示等）。 您需要 totalcount 和 lastmedian，以及 newvalue。
{
totalcount++;
newmedian=lastmedian+(newvalue>lastmedian?1:-1)*(lastmedian==0?newvalue: lastmedian/totalcount*2);
}
对于像页面显示时间这样的事情，可以产生相当精确的结果。
规则：输入流需要在页面显示时间的顺序上平稳，数量大（>30等），并且具有非零中位数。
例如： 页面加载时间，800个项目，10毫秒...3000毫秒，平均90毫秒，真实中位数：11毫秒
经过30次输入后，中位数误差通常<=20%（9ms..12ms），并且越来越少。 经过800次输入后，误差为+-2%。
另一个思考者提出了类似的解决方案：Median Filter Super efficient implementation