我目前正在研究一种在 C 中实现滚动中值滤波器（类似于滚动均值滤波器）的算法。根据我对文献的搜索，似乎有两种相当有效的方法可以做到这一点。第一个是对初始值窗口进行排序，然后执行二进制搜索以插入新值并在每次迭代时删除现有值。
第二个（来自 Hardle 和 Steiger，1995，JRSS-C，算法 296）构建了一个双端堆结构，一端是 maxheap，另一端是 minheap，中间是中间值。这产生了一种线性时间算法，而不是 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
我看过 Rsrc/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 建议的顶级编码器链接中实现了这两个想法（匹配编辑：向下滚动到 FloatingMedian）。
第一个想法将数据划分为两个数据结构（堆、多集等），每个插入/删除 O(ln N) 不允许在没有大成本的情况下动态更改分位数。即我们可以有一个滚动的中位数，或者一个滚动的 75%，但不能同时有两者。
第二个想法使用 O(ln N) 的段树进行插入/删除/查询，但更灵活。最好的“N”是数据范围的大小。因此，如果您的滚动中位数具有一百万个项目的窗口，但您的数据与 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（搜索“拆分和连接”）。
我已经包装了树，以便在支持 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 实现。这个问题还有更多细节：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
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
滚动中位数可以通过维护两个数字分区来找到。
维护分区使用最小堆和最大堆。
最大堆将包含小于等于中位数的数字。
最小堆将包含大于等于中位数的数字。
平衡约束： 如果元素总数是偶数，那么两个堆应该有相等的元素。
如果元素总数为奇数，则 Max Heap 将比 Min Heap 多一个元素。
中值元素：如果两个分区具有相同数量的元素，则中值将是第一个分区的最大元素和第二个分区的最小元素之和的一半。
否则中位数将是第一个分区的最大元素。
算法- 1-取两个堆（1个最小堆和1个最大堆） Max Heap 将包含前半数的元素 最小堆将包含后半数的元素 2-将流中的新数字与最大堆顶部进行比较， 如果它更小或相等，则将该数字添加到最大堆中。 否则在最小堆中添加数字。 3- 如果 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 个元素和一个计数的数组，然后清除所有这些值。对于流中的每个值，递增相应的 bin 和计数。在流结束后找到包含计数/2 最高值的 bin - 通过添加从 0 开始的连续 bin 轻松完成。使用相同的方法可以找到任意排名顺序的值。（如果需要在运行期间检测箱饱和并将存储箱的大小“升级”为更大的类型，则会有一个小麻烦。）
这种特殊情况看似人为，但在实践中却很常见。如果实数在一个范围内并且已知“足够好”的精度水平，它也可以用作实数的近似值。这几乎适用于对一组“现实世界”对象的任何一组测量。例如，一群人的身高或体重。集数不够大？它同样适用于地球上所有（个体）细菌的长度或重量——假设有人可以提供数据！
看起来我误读了原文——它似乎想要一个滑动窗口中位数，而不是一个很长的流的中位数。这种方法仍然适用。为初始窗口加载前 N 个流值，然后为第 N+1 个流值增加相应的 bin，同时减少与第 0 个流值对应的 bin。在这种情况下，有必要保留最后 N 个值以允许递减，这可以通过循环寻址大小为 N 的数组来有效地完成。由于中位数的位置只能改变 -2,-1,0,1 ,2 在滑动窗口的每一步上，没有必要将所有 bin 相加到每一步的中值，只需根据修改了哪些边 bin 来调整“中值指针”。例如，如果新值和被删除的值都低于当前中值，则它不会改变（偏移量 = 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 以及 newvalue。
{
totalcount++;
newmedian=lastmedian+(newvalue>lastmedian?1:-1)*(lastmedian==0?newvalue: lastmedian/totalcount*2);
}
为诸如 page_display_time 之类的内容生成非常精确的结果。
规则：输入流需要按页面显示时间顺序平滑，数量大（>30 等），并且中值不为零。
示例：页面加载时间，800 项，10ms...3000ms，平均 90ms，实际中位数：11ms
30 次输入后，中位数误差一般 <=20% (9ms..12ms)，并且越来越小。输入 800 次后，误差为 +-2%。
另一个有类似解决方案的思想家在这里：中值滤波器超高效实现
基于@mathog 的想法，这是一个 C# 实现，用于在具有已知值范围的字节数组上运行中位数。可以扩展到其他整数类型。
/// <summary>
/// Median estimation by histogram, avoids multiple sorting operations for a running median
/// </summary>
public class MedianEstimator
{
private readonly int m_size2;
private readonly byte[] m_counts;
/// <summary>
/// Estimated median, available right after calling <see cref="Init"/> or <see cref="Update"/>.
/// </summary>
public byte Median { get; private set; }
/// <summary>
/// Ctor
/// </summary>
/// <param name="size">Median size in samples</param>
/// <param name="maxValue">Maximum expected value in input data</param>
public MedianEstimator(
int size,
byte maxValue)
{
m_size2 = size / 2;
m_counts = new byte[maxValue + 1];
}
/// <summary>
/// Initializes the internal histogram with the passed sample values
/// </summary>
/// <param name="values">Array of values, usually the start of the array for a running median</param>
public void Init(byte[] values)
{
for (var i = 0; i < values.Length; i++)
m_counts[values[i]]++;
UpdateMedian();
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
private void UpdateMedian()
{
// The median is the first value up to which counts add to size / 2
var sum = 0;
Median = 0;
for (var i = 0; i < m_counts.Length; i++)
{
sum += m_counts[i];
Median = (byte) i;
if (sum > m_size2) break;
}
}
/// <summary>
/// Updates the median estimation by removing <paramref name="last"/> and adding <paramref name="next"/>. These
/// values must be updated as the running median is applied. If the median length is <i>N</i>, at the sample
/// <i>i</i>, <paramref name="last"/> is sample at index <i>i</i>-<i>N</i>/2 and <paramref name="next"/> is sample
/// at index <i>i</i>+<i>N</i>/2+1.
/// </summary>
/// <param name="last">Sample at the start of the moving window that is to be removed</param>
/// <param name="next">Sample at the end of the moving window + 1 that is to be added</param>
public void Update(byte last, byte next)
{
m_counts[last]--;
m_counts[next]++;
// The conditions below do not change median value so there is no need to update it
if (last == next ||
last < Median && next < Median || // both below median
last > Median && next > Median) // both above median
return;
UpdateMedian();
}
针对运行中位数进行测试，有时间：
private void TestMedianEstimator()
{
var r = new Random();
const int SIZE = 15;
const byte MAX_VAL = 80;
var values = new byte[100000];
for (var i = 0; i < values.Length; i++)
values[i] = (byte) (MAX_VAL * r.NextDouble());
var timer = Stopwatch.StartNew();
// Running median
var window = new byte[2 * SIZE + 1];
var medians = new byte[values.Length];
for (var i = SIZE; i < values.Length - SIZE - 1; i++)
{
for (int j = i - SIZE, k = 0; j <= i + SIZE; j++, k++)
window[k] = values[j];
Array.Sort(window);
medians[i] = window[SIZE];
}
timer.Stop();
var elapsed1 = timer.Elapsed;
timer.Restart();
var me = new MedianEstimator(2 * SIZE + 1, MAX_VAL);
me.Init(values.Slice(0, 2 * SIZE + 1));
var meMedians = new byte[values.Length];
for (var i = SIZE; i < values.Length - SIZE - 1; i++)
{
meMedians[i] = me.Median;
me.Update(values[i - SIZE], values[i + SIZE + 1]);
}
timer.Stop();
var elapsed2 = timer.Elapsed;
WriteLineToLog($"{elapsed1.TotalMilliseconds / elapsed2.TotalMilliseconds:0.00}");
var diff = 0;
for (var i = 0; i < meMedians.Length; i++)
diff += Math.Abs(meMedians[i] - medians[i]);
WriteLineToLog($"Diff: {diff}");
}
这是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)，然后将它们相加。然后这成为新的平均值。
编辑：不是用户要求的，也不是统计上有效的，但对于很多用途来说已经足够了。
我会把它留在这里（尽管有反对票）以供搜索！