【问题标题】:how to calculate running median efficiently如何有效地计算运行中位数
【发布时间】:2016-10-06 21:18:31
【问题描述】:

我借用了一些代码,试图实现一个函数来计算大量数据的运行中位数。当前的对我来说太慢了(棘手的部分是我需要从运行框中排除所有零)。下面是代码:

from itertools import islice
from collections import deque
from bisect import bisect_left,insort

def median(s):
    sp = [nz for nz in s if nz!=0]
    print sp
    Mnow = len(sp)
    if Mnow == 0:
        return 0
    else:
        return np.median(sp)

def RunningMedian(seq, M):
    seq = iter(seq)
    s = []

    # Set up list s (to be sorted) and load deque with first window of seq
    s = [item for item in islice(seq,M)]
    d = deque(s)

    # Sort it in increasing order and extract the median ("center" of the sorted window)
    s.sort()
    medians = [median(s)]
    for item in seq:
        old = d.popleft()          # pop oldest from left
        d.append(item)             # push newest in from right
        del s[bisect_left(s, old)] # locate insertion point and then remove old 
        insort(s, item)            # insert newest such that new sort is not required        
        medians.append(median(s))
    return medians

效果很好,唯一的缺点就是太慢了。任何人都可以帮助我改进代码以提高效率?

在我探索了所有可能性之后,下面的简单代码可以比较高效地计算:

def RunningMedian(x,N):
    idx = np.arange(N) + np.arange(len(x)-N+1)[:,None]
    b = [row[row>0] for row in x[idx]]
    return np.array(map(np.median,b))
    #return np.array([np.median(c) for c in b])  # This also works

【问题讨论】:

  • 如果您将此问题移至 Code Review Stack Exchange,您可以获得更多答案。
  • Running or sliding median, mean and standard deviation 不适合您,还是太慢了,无法满足您的需求?
  • 看看this。这是很多代码,用 Python 编写可能不会那么快。但是实现的速度差不多,它使用与 pandas 进行滚动中位数计算相同的方法,但为了速度是用 Cython 编写的。
  • 为什么不使用 pandas.rolling_* 方法?

标签: python numpy slice median


【解决方案1】:

一种方法如下:

def RunningMedian(x,N):
    idx = np.arange(N) + np.arange(len(x)-N+1)[:,None]
    b = [row[row>0] for row in x[idx]]
    return np.array(map(np.median,b))
    #return np.array([np.median(c) for c in b])  # This also works

我找到了一个更快的(快几万倍),复制如下:

from collections import deque
from bisect import insort, bisect_left
from itertools import islice
def running_median_insort(seq, window_size):
    """Contributed by Peter Otten"""
    seq = iter(seq)
    d = deque()
    s = []
    result = []
    for item in islice(seq, window_size):
        d.append(item)
        insort(s, item)
        result.append(s[len(d)//2])
    m = window_size // 2
    for item in seq:
        old = d.popleft()
        d.append(item)
        del s[bisect_left(s, old)]
        insort(s, item)
        result.append(s[m])
    return result

看看链接:running_median

【讨论】:

  • "A much faster one" 当窗口长度是偶数时给出不正确的中位数。正确的做法是在一对中心元素之间取均值。
  • 您可以检查一些简单的情况,例如窗口大小为 3 的阶跃函数 - 这两个函数在阶跃函数的边界上给出了不正确的结果。而且两者都不保留相位。
【解决方案2】:

您可以在处理数据样本时使用两个heaps 来维护数据样本的下半部分和上半部分。算法是这样的:对于每个值,将其放入适当的堆并“重新平衡”堆,以使它们的大小相差不超过 1。然后,要获得中位数,只需从较大的元素中拉出第一个元素堆,或者如果两个堆的大小相等,则取两个堆的第一个元素的平均值。该解决方案的时间复杂度为O(n log(n))

from heapq import heappush, heappop


class RunningMedian:
    def __init__(self):
        self.lowers, self.highers = [], []

    def add_number(self, number):
        if not self.highers or number > self.highers[0]:
            heappush(self.highers, number)
        else:
            heappush(self.lowers, -number)  # for lowers we need a max heap
        self.rebalance()

    def rebalance(self):
        if len(self.lowers) - len(self.highers) > 1:
            heappush(self.highers, -heappop(self.lowers))
        elif len(self.highers) - len(self.lowers) > 1:
            heappush(self.lowers, -heappop(self.highers))

    def get_median(self):
        if len(self.lowers) == len(self.highers):
            return (-self.lowers[0] + self.highers[0])/2
        elif len(self.lowers) > len(self.highers):
            return -self.lowers[0]
        else:
            return self.highers[0]

演示:

>>> running_median = RunningMedian()
>>> for n in (12, 4, 5, 3, 8, 7):
...     running_median.add_number(n)
...     print(running_median.get_median())
... 
12
8.0
5
4.5
5
6.0

【讨论】:

  • 这是怎么回事?堆无限增长。通常,您在样本上有一个窗口来计算运行中位数。因此,您既需要在窗口中订购样本,又需要能够按照它们被引入的顺序弹出样本。一个有效的算法应该在添加/删除元素时重用排序——你的实现会这样做——同时仍然保持窗口大小——你的实现不会。
猜你喜欢
  • 2012-06-11
  • 1970-01-01
  • 2010-11-13
  • 2016-03-22
  • 1970-01-01
  • 1970-01-01
  • 2015-05-19
  • 1970-01-01
  • 2014-06-14
相关资源
最近更新 更多