how to calculate running median efficiently

问题:

I borrowed some code trying to implement a function to calculate the running median for a ton of data. The current one is too slow for me (The tricky part is that I need to exclude all zeros from the running box). Below is the code:

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

It works well, the only drawback is that it is too slow. Any one could help me to improve the code to be more efficient? Thanks.

After I explored all the possibilities, the following simple code can calculate comparably efficiently:

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

Thank @Divakar.

回答1:

One approach is below:

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

I found a much faster one (tens of thousand times faster), copied as below:

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

Take a look at the link: running_median

回答2:

One way is to use two heaps to maintain the lower and higher halves of the data sample as you process it. The algorithm: for each value, put it to an appropriate heap and 'rebalance' the heaps so that their size is not different by more than 1. Then to get the median just pull out the first element from the bigger heap or take an average of the first elements of the two heaps. This solutions has O(n log(n)) time complexity.

from heapq import heappush, heappop


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


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


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


lowers, highers = [], []

for i in [12, 4, 5, 3, 8, 7]:
    add_number(i, lowers, highers)
    rebalance(lowers, highers)
    print(get_median(lowers, highers))