那些阅读过我之前的问题的人都知道我在理解和实现快速排序和快速选择以及其他一些基本算法方面所做的工作。
快速选择用于计算未排序列表中的第k个最小元素，这个概念也可以用于查找未排序列表中的中位数。
这一次，我需要帮助设计一种有效的技术来计算运行中位数，因为 quickselect 不是一个好的选择，因为它需要在每次列表更改时重新计算。因为 quickselect 每次都必须重新启动，它不能利用之前完成的计算，所以我正在寻找一种不同的算法，它相似（可能）但在运行中位数方面更有效。
【问题讨论】：
#include<cstdio>
#include<iostream>
#include<queue>
#include <vector>
#include <functional>
typedef priority_queue<unsigned int> type_H_low;
typedef priority_queue<unsigned int, std::vector<unsigned int>, std::greater<unsigned int> > type_H_high;
size_t signum(int left, int right) {
if (left == right){
return 0;
}
return (left < right)?-1:1;
}
void get_median( unsigned int x_i, unsigned int &m, type_H_low *l, type_H_high *r) {
switch (signum( l->size(), r->size() )) {
case 1: // There are more elements in left (max) heap
if (x_i < m) {
r->push(l->top());
l->pop();
l->push(x_i);
} else {
r->push(x_i);
}
break;
case 0: // The left and right heaps contain same number of elements
if (x_i < m){
l->push(x_i);
} else {
r->push(x_i);
}
break;
case -1: // There are more elements in right (min) heap
if (x_i < m){
l->push(x_i);
} else {
l->push(r->top());
r->pop();
r->push(x_i);
}
break;
}
if (l->size() == r->size()){
m = l->top();
} else if (l->size() > r->size()){
m = l->top();
} else {
m = r->top();
}
return;
}
void print_median(vector<unsigned int> v) {
unsigned int median = 0;
long int sum = 0;
type_H_low H_low;
type_H_high H_high;
for (vector<unsigned int>::iterator x_i = v.begin(); x_i != v.end(); x_i++) {
get_median(*x_i, median, &H_low, &H_high);
std::cout << median << std::endl;
}
}
【讨论】：
滚动窗口中位数算法：
median 是一个排序数组，您可以从中获取中间值。
简单的滚动实现是一个队列（dqueue）和一个 sorted_array（任何实现，二叉树，skiparray）。
d_queue 是一个数组，您可以在其中推到尾部并从数组的前面移位（弹出）。
sorted_array 是一个数组，您可以在其中按顺序插入使用二进制搜索找到的位置。
我使用队列（先进先出数组）来跟踪添加值的顺序，以便在队列长于所需大小时知道要从中间数组中删除哪些项目。要按日期时间或某个运行索引删除元素，可以添加另一个队列并检查第一个元素是否太旧，然后决定是否从两个队列中删除第一个值。
为了有效地计算中位数，我使用了排序数组技术。当您将新项目插入其排序位置时，数组始终是排序的。
插入：
删除：
要有中位数：
【讨论】：
这是一个 C++ 平衡树结构，它提供了按排序列表中的索引进行查询的能力。由于它按排序顺序维护所有值，因此它不如双堆方法高效，但它提供了一些额外的灵活性。例如，它还可以为您提供跑步四分位数。
template <typename T>
class Node
{
public:
T key;
Node* left;
Node* right;
size_t size;
Node(T k) : key(k)
{
isolate();
}
~Node()
{
delete(left);
delete(right);
}
void isolate()
{
left = NULL;
right = NULL;
size = 1;
}
void recount()
{
size = 1 + (left ? left->size : 0) + (right ? right->size : 0);
}
Node<T>* rotateLeft()
{
Node<T>* c = right;
Node<T>* gc = right->left;
right = gc;
c->left = this;
recount();
c->recount();
return c;
}
Node<T>* rotateRight()
{
Node<T>* c = left;
Node<T>* gc = left->right;
left = gc;
c->right = this;
recount();
c->recount();
return c;
}
Node<T>* balance()
{
size_t lcount = left ? left->size : 0;
size_t rcount = right ? right->size : 0;
if((lcount + 1) * 2 < (rcount + 1))
{
size_t lcount2 = right->left ? right->left->size : 0;
size_t rcount2 = right->right ? right->right->size : 0;
if(lcount2 > rcount2)
right = right->rotateRight();
return rotateLeft();
}
else if((rcount + 1) * 2 <= (lcount + 1))
{
size_t lcount2 = left->left ? left->left->size : 0;
size_t rcount2 = left->right ? left->right->size : 0;
if(lcount2 < rcount2)
left = left->rotateLeft();
return rotateRight();
}
else
{
recount();
return this;
}
}
Node<T>* insert(Node<T>* newNode)
{
if(newNode->key < key)
{
if(left)
left = left->insert(newNode);
else
left = newNode;
}
else
{
if(right)
right = right->insert(newNode);
else
right = newNode;
}
return balance();
}
Node<T>* get(size_t index)
{
size_t lcount = left ? left->size : 0;
if(index < lcount)
return left->get(index);
else if(index > lcount)
return right ? right->get(index - lcount - 1) : NULL;
else
return this;
}
Node<T>* find(T k, size_t start, size_t* outIndex)
{
if(k < key)
return left ? left->find(k, start, outIndex) : NULL;
else if(k > key)
return right ? right->find(k, left ? start + left->size + 1 : start + 1, outIndex) : NULL;
else
{
if(outIndex)
*outIndex = start + (left ? left->size : 0);
return this;
}
}
Node<T>* remove_by_index(size_t index, Node<T>** outNode)
{
size_t lcount = left ? left->size : 0;
if(index < lcount)
left = left->remove_by_index(index, outNode);
else if(index > lcount)
right = right->remove_by_index(index - lcount - 1, outNode);
else
{
*outNode = this;
size_t rcount = right ? right->size : 0;
if(lcount < rcount)
return left ? right->insert(left) : right;
else
return right ? left->insert(right) : left;
}
return balance();
}
Node<T>* remove_by_value(T k, Node<T>** outNode)
{
if(k < key)
{
if(!left)
throw "not found";
left = left->remove_by_value(k, outNode);
}
else if(k > key)
{
if(!right)
throw "not found";
right = right->remove_by_value(k, outNode);
}
else
{
*outNode = this;
size_t lcount = left ? left->size : 0;
size_t rcount = right ? right->size : 0;
if(lcount < rcount)
return left ? right->insert(left) : right;
else
return right ? left->insert(right) : left;
}
return balance();
}
};
template <typename T>
class MyReasonablyEfficientRunningSortedIndexedCollection
{
private:
Node<T>* root;
Node<T>* spare;
public:
MyReasonablyEfficientRunningSortedIndexedCollection() : root(NULL), spare(NULL)
{
}
~MyReasonablyEfficientRunningSortedIndexedCollection()
{
delete(root);
delete(spare);
}
void insert(T key)
{
if(spare)
spare->key = key;
else
spare = new Node<T>(key);
if(root)
root = root->insert(spare);
else
root = spare;
spare = NULL;
}
void drop_by_index(size_t index)
{
if(!root || index >= root->size)
throw "out of range";
delete(spare);
root = root->remove_by_index(index, &spare);
spare->isolate();
}
void drop_by_value(T key)
{
if(!root)
throw "out of range";
delete(spare);
root = root->remove_by_value(key, &spare);
spare->isolate();
}
T get(size_t index)
{
if(!root || index >= root->size)
throw "out of range";
return root->get(index)->key;
}
size_t find(T key)
{
size_t outIndex;
Node<T>* node = root ? root->find(key, 0, &outIndex) : NULL;
if(node)
return outIndex;
else
throw "not found";
}
size_t size()
{
return root ? root->size : 0;
}
};
【讨论】：
Jeff McClintock 运行中位数估计。只需要保留两个值。 此示例迭代一组采样值（CPU 消耗）。似乎相对较快地收敛（大约 100 个样本）到中位数的估计值。 这个想法是在每次迭代中，中位数英寸以恒定速率朝向输入信号。该比率取决于您估计中位数的大小。我使用平均值作为中位数大小的估计，以确定中位数每个增量的大小。如果您需要精确到大约 1% 的中位数，请使用 0.01 * 平均值的步长。
float median = 0.0f;
float average = 0.0f;
// for each sample
{
average += ( sample - average ) * 0.1f; // rough running average.
median += _copysign( average * 0.01, sample - median );
}
【讨论】：
average * 0.01 和采样率频率为 3) 的抖动，在短脉冲上偏转（中值原本没有，因此作为胡椒和噪声滤波器很受欢迎）
streaming median 是使用两个堆计算的。所有小于或等于当前中位数的数都在左堆中，左堆的排列使得最大数在堆的根。所有大于或等于当前中位数的数字都在右堆中，它的排列使得最小数字在堆的根。请注意，等于当前中位数的数字可以在任一堆中。两个堆中的数字计数相差永远不会超过 1。
当进程开始时，两个堆最初是空的。输入序列中的第一个数字被添加到其中一个堆中，无论哪个都没有关系，并作为第一个流媒体中位数返回。然后将输入序列中的第二个数字添加到另一个堆，如果右堆的根小于左堆的根，则交换两个堆，并将两个数字的平均值作为第二个流返回中位数。
然后主算法开始。输入序列中的每个后续数字都与当前中位数进行比较，如果小于当前中位数则添加到左堆，如果大于当前中位数则添加到右堆；如果输入数等于当前中位数，则将其添加到具有较小计数的堆中，或者如果它们具有相同的计数，则任意添加到任一堆中。如果这导致两个堆的计数相差超过 1，则删除较大堆的根并插入较小的堆。然后当前中位数被计算为较大堆的根，如果它们的计数不同，或者两个堆的根的平均值，如果它们大小相同。
Scheme 和 Python 中的代码可在我的blog 获得。
【讨论】：
一种解决方案是维护一个order statistic tree，依次插入序列中的每个元素，然后计算树中元素的中位数。
每次插入需要 O(lg n) 时间，每个中位数需要 O(lg n) 时间，总共需要 O(n em> lg n) 时间，加上 O(n) 空间。
【讨论】：