how to calculate running median efficiently

Question

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?

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

Possible duplicate of [Rolling median algorithm in C](http://stackoverflow.com/questions/1309263/rolling-median-algorithm-in-c) — Josephine Moeller, Jun 07 '16 at 05:57
You can attract more answers if you move this question to Code Review Stack Exchange. — Selcuk, Jun 07 '16 at 06:44
Didn't [Running or sliding median, mean and standard deviation](http://stackoverflow.com/a/33585850/3293881) work for you or was it too slow for your needs? — Divakar, Jun 07 '16 at 06:55
@Divakar Thank you very much. It works, but still too slow for a huge dataset. — Huanian Zhang, Jun 07 '16 at 06:58
@Divakar And also I want to exclude the zeros from the running box. I did not how to implement that from your answer. Do you have any ideas? Thank you. — Huanian Zhang, Jun 07 '16 at 07:11
What is *too slow* for you? How big is your data and how long do you expect it to take? Median of large datasets is a quite expensive calculation... a rolling median is even more. — Imanol Luengo, Jun 07 '16 at 09:24
Take a look at [this](http://code.activestate.com/recipes/576930-efficient-running-median-using-an-indexable-skipli/). It's a lot of code, and being written in Python may not be all that fast. But the implementation is about as fast as it gets, and it uses the same approach that pandas does for its rolling median calculations, but written in Cython for speed. — Jaime, Jun 07 '16 at 12:34
@ImanolLuengo The array has 300 elements, I want the window size to be 25, and there are 20 million such arrays to be analyzed. — Huanian Zhang, Jun 07 '16 at 22:47
@Jaime I see some similar codes, but they do not have the requirement that the running box should exclude all zeros. This is the tricky part. — Huanian Zhang, Jun 07 '16 at 22:48
@Apero Again the pandas can not deal with non-zero running box. — Huanian Zhang, Jun 07 '16 at 22:49

Huanian Zhang · Accepted Answer · 2017-08-21T22:37:52.180

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

"A much faster one" gives incorrect median when the length of the window is even. Correct one is to take mean between a pair of center elements. — Maksym Ganenko, Aug 07 '19 at 13:13
You may check some simple case like step function with window size 3 - both functions give incorrect results on boundary of the step function. And both doesn't preserve phase. — Maksym Ganenko, Aug 07 '19 at 13:24

Eugene Yarmash · Answer 2 · 2019-08-20T15:32:06.463

You can use two heaps to maintain the lower and higher halves of the data sample as you process it. The algorithm goes like this: 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 if they are of equal size. This solution has O(n log(n)) time complexity.

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]

Demo:

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

What is running about this? The heap grows without bounds. Normally you have a window over the samples to compute the running median. Therefore you need both to order the samples in the window and to be able to pop samples in the order that they have been introduced. An efficient algorithm should reuse the ordering when adding/removing elements - which your implementation does - while still keeping window size - which your implementation does not. — blueFast, Sep 23 '20 at 09:20

how to calculate running median efficiently

2 Answers2

Linked

Related