C++高效计算运行中位数

27

如果您已经阅读过我的以前的问题,您就会知道我在理解和实现快速排序、快速选择和其他基本算法方面的工作。

快速选择用于计算未排序列表中第k个最小元素,这个概念也可以用于查找未排序列表中的中位数。

这一次,我需要帮助设计一种高效的技术来计算滑动中位数,因为快速选择不是一个好的选择,它需要每次列表更改时重新计算。因为快速选择必须每次重新开始,所以无法利用以前的计算结果,因此我正在寻找一种类似但更有效的算法来解决滑动中位数领域的问题。


可以使用快速排序算法中的分区来在线性时间内完成,但最坏情况下需要n^2的时间。在您的集合中选择随机点作为枢轴,并移动其他元素,使小于枢轴的元素位于左侧,大于或等于枢轴的元素位于右侧。如果枢轴在中间,则它是中位数,否则转到具有中位数(较大尺寸块)的块。重复。另一个保证线性时间的算法是中位数中的中位数,CLRS中有描述,我相信维基百科上也有相关内容。请查阅。 - Adrian
6个回答

51

使用两个堆来计算流媒体中位数。所有小于或等于当前中位数的数字都在左堆中,该堆被排列使得最大数字位于堆的根部。所有大于或等于当前中位数的数字都在右堆中,该堆被排列使得最小数字位于堆的根部。请注意,等于当前中位数的数字可以在任一堆中。两个堆中数字的数量永远不会相差超过1。

当进程开始时,两个堆最初为空。将输入序列中的第一个数字添加到其中一个堆中,无论哪个都可以,并返回作为第一个流媒体中位数。然后将输入序列中的第二个数字添加到另一个堆中,如果右堆的根小于左堆的根,则交换两个堆,并返回两个数字的平均值作为第二个流媒体中位数。

然后主算法开始。将输入序列中的每个后续数字与当前中位数进行比较,如果小于当前中位数,则将其添加到左堆中;如果大于当前中位数,则将其添加到右堆中;如果输入数字等于当前中位数,则将其添加到计数较小的堆中,或者如果它们的计数相同,则随意添加到任一堆中。如果这导致两个堆的计数相差超过1,则从较大的堆中删除根并将其插入较小的堆中。然后计算当前中位数,如果它们的计数不同,则为较大堆的根,如果它们的大小相同,则为两个堆的根的平均值。

在我的blog上提供了Scheme和Python代码。


1
有没有用C++实现的代码?顺便感谢您的回复,非常感谢。我喜欢详细的解释。 - Edge
另外,你的想法只适用于排序列表还是未排序列表也可以? - Edge
@Andrew,你看过boost累加器吗? - Chisholm
我不了解boost函数。 - Edge
4
添加元素到系统中是可以的。但是如何删除旧元素呢? - Notinlist
显示剩余5条评论

19

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

enter image description here


2
虽然这个解决方案非常高效,但请注意以下注意事项:1)收敛速度取决于信号幅度(比较具有不同偏移和幅度的阶跃响应),因此不会收敛到接近零的信号!2)对于接近恒定输入信号,该估计引入了振荡,其振幅为“平均值* 0.01”,频率为采样率3)在短脉冲上偏转(原始中位数不会这样做,因此作为胡椒和噪声滤波器而受欢迎) - orzechow
使用滚动标准差而不是滚动平均值来缩放步增量可能是针对具有对数正态分布的数据(这种数据占许多/大多数自然发生的过程)的有趣方向。使用变异性度量而不是平均值将解决Orzechow提出的抖动问题。不幸的是,方差不具有稳健性,因此大的瞬时异常值可能会重新引入抖动。然后,叠加滤波器可能会产生增加有效窗口宽度的效果,这很有趣。 - sircolinton

6
一种解决方案是维护一个有序统计树,将序列中的每个元素依次插入,然后计算树中元素的中位数。
每次插入需要 O(lg n) 的时间,每次计算中位数也需要 O(lg n) 的时间,总共需要 O(n lg n) 的时间,加上 O(n) 的空间。

这种树类型对于这个目的来说好吗?我以前没听说过它。 - Edge
2
有序统计树允许在对数时间内进行索引,即获取序列的第k小元素。 - Fred Foo
这个能在未排序的列表中工作吗? - Edge
其实算了,对我来说看起来很简单。还有其他参考资料可以更深入地讨论顺序统计树的工作原理吗? - Edge
@Andrew:https://zh.wikipedia.org/wiki/%E7%AE%97%E6%B3%95%E5%AF%BC%E8%AE%BA - Fred Foo
显示剩余2条评论

1

滚动窗口中位数算法:

中位数是一个排序数组,你从中取出中间值。

简单的滚动实现可以使用队列(双端队列)和一个排序数组(任何实现,如二叉树、跳表)。

双端队列是一个数组,你可以在尾部插入(push)和在头部移除(shift/pop)数组元素。

排序数组是一个数组,你可以通过二分查找找到要插入的位置并按顺序插入。

我使用了队列(先进先出的数组)来跟踪添加值的顺序,以便知道哪些项目从中位数数组中删除,当队列长度超过所需大小时。为了按日期时间或某个运行索引掉落元素,可以添加另一个队列,并检查第一个元素是否太旧,并决定是否从两个队列中同时移除第一个值。

为了有效地计算中位数,我使用了一种排序数组技术。当你将新项插入其排序位置时,数组始终保持排序状态。

  1. 插入:

    • 在有序数组中的有序位置插入值,
    • 并将一个值推入队列。
  2. 删除:

    • 如果d_queue的第一个元素超出了窗口,或者如果在另一个具有索引的队列中,索引太旧,则:
      • 从d_queue(s)中删除第一个项目,
      • 并在排序数组中进行二进制搜索并删除它。
  3. 获取中位数:

    • 使用排序数组中间的值。
    • 如果排序数组长度为偶数,则使用中间项。
    • 如果排序数组长度为奇数,则使用两个中间项的平均值。

1
这是一个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;
    }
};

0
#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;
    }
}

2
你好,欢迎来到SO。你的答案只包含代码。如果你能加一些注释来解释它是如何工作的,那就更好了。你能否请编辑你的答案并添加注释?谢谢! - Fabio says Reinstate Monica

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接