C ++有效地计算运行中位数[重复]

Posted

技术标签:

【中文标题】C ++有效地计算运行中位数[重复]【英文标题】：C++ Efficiently Calculating a Running Median [duplicate] 【发布时间】：2012-06-11 10:39:01 【问题描述】：

那些阅读过我之前的问题的人都知道我在理解和实现快速排序和快速选择以及其他一些基本算法方面所做的工作。

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

这一次，我需要帮助设计一种有效的技术来计算运行中位数，因为 quickselect 不是一个好的选择，因为它需要在每次列表更改时重新计算。因为 quickselect 每次都必须重新启动，它不能利用之前完成的计算，所以我正在寻找一种不同的算法，它相似（可能）但在运行中位数方面更有效。

【问题讨论】：

这可以通过使用快速排序算法中的分区在线性时间内完成，但时间最差为 n^2。选择集合中的随机点作为枢轴并移动其他元素，使小于枢轴的元素在左侧，大于或等于的元素在右侧。如果枢轴在中间，则为中值，如果不是，则转到具有中值的块（较大的块）。重复。另一种保证线性时间的算法，它是 CLRS 中描述的中位数的中位数，我也相信在***上。查一下。 【参考方案1】： 
#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;

【讨论】：

您好，欢迎来到 SO。您的答案仅包含代码。如果您还可以添加一些评论来解释它的作用和方式，那就更好了。你能请edit你的答案并添加吗？谢谢！【参考方案2】：

滚动窗口中位数算法：

median 是一个排序数组，您可以从中获取中间值。

简单的滚动实现是一个队列（dqueue）和一个 sorted_array（任何实现，二叉树，skiparray）。

d_queue 是一个数组，您可以在其中推到尾部并从数组的前面移位（弹出）。

sorted_array 是一个数组，您可以在其中按顺序插入使用二进制搜索找到的位置。

我使用队列（先进先出数组）来跟踪添加值的顺序，以便在队列长于所需大小时知道要从中间数组中删除哪些项目。要按日期时间或某个运行索引删除元素，可以添加另一个队列并检查第一个元素是否太旧，然后决定是否从两个队列中删除第一个值。

为了有效地计算中位数，我使用了排序数组技术。当您将新项目插入其排序位置时，数组始终是排序的。

    插入：

    在 sorted_array 中的有序位置插入， 并将值推送到队列中。

    删除：

    如果 d_queue 第一个元素不在窗口中，或者如果在另一个队列中可以有索引，则索引太旧，则： 从 d_queue(s) 中删除第一项， 并在排序后的数组中对其进行二进制搜索并将其删除。

    要有中位数：

    使用 sorted_array 中间的值。 如果 sorted_array 长度甚至使用中间的项目。 如果 sorted_array 长度为奇数，则在中间使用两项的平均值。

【讨论】：

【参考方案3】：

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

【讨论】：

【参考方案4】：

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

【讨论】：

虽然此解决方案非常有效，但请注意以下警告：1) 收敛速度取决于信号幅度（比较具有不同偏移和幅度的阶跃响应），因此不会针对接近零的信号收敛！ 2) 对于接近恒定的输入信号，该估计引入了幅度为average * 0.01 和采样率频率为 3) 的抖动，在短脉冲上偏转（中值原本没有，因此作为胡椒和噪声滤波器很受欢迎） 使用滚动标准差而不是滚动平均值来缩放步长增量可能是对数正态分布数据的一个有趣方向（它解释了许多/最自然产生的过程）。使用可变性度量而不是平均值将解决由 orzechow 引起的抖动问题。不幸的是，方差并不稳健，因此大的瞬态异常值可能会重新引入抖动。然后堆叠过滤器可能会很有趣，我希望它具有增加有效窗口宽度的效果。【参考方案5】：

streaming median 是使用两个堆计算的。所有小于或等于当前中位数的数都在左堆中，左堆的排列使得最大数在堆的根。所有大于或等于当前中位数的数字都在右堆中，它的排列使得最小数字在堆的根。请注意，等于当前中位数的数字可以在任一堆中。两个堆中的数字计数相差永远不会超过 1。

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

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

Scheme 和 Python 中的代码可在我的blog 获得。

【讨论】：

有C++实现的代码吗？顺便谢谢您的回复，不胜感激。我喜欢详细的解释。 另外，你的想法是只适用于排序列表还是未排序列表？ @Andrew，你看过升压蓄电池吗？ 我对 boost 函数一无所知。 向系统添加元素是可以的。但是如何删除旧元素呢？【参考方案6】：

一种解决方案是维护一个order statistic tree，依次插入序列中的每个元素，然后计算树中元素的中位数。

每次插入需要 O(lg n) 时间，每个中位数需要 O(lg n) 时间，总共需要 O(n em> lg n) 时间，加上 O(n) 空间。

【讨论】：

这种树型适合这个目的吗？我以前没听说过。 顺序统计树允许索引，即在对数时间内获取序列的第 k 个最小元素。 这是否适用于未排序的列表？ 是的，您只需将项目插入树中，它们就会被排序。顺序统计树是增强的 BST。 我注意到您链接的页面包含伪代码，而不是实际代码。我了解基本的 BST，但仅此而已，并且在伪代码之后实现该树需要一段时间。

以上是关于C ++有效地计算运行中位数[重复]的主要内容，如果未能解决你的问题，请参考以下文章

运行或滑动中位数、平均值和标准差

如何计算 PrestoSQL 中的中位数？

在列表中有效地重复data.table，从循环中的另一个data.table顺序替换具有相同名称的列

计算 TB 数据集中分位数的高效算法

如何在python中有效地计算（稀疏）位矩阵的矩阵乘积

在 python 中有效地处理一个大的 .txt 文件