# Algorithms to compute the running median?

On smaller window sizes, n log n sorting might work. Are there any better algorithms to achieve this?

• I think this is the first candidate to be moved to Stack Overflow. – user88 Jul 19 '10 at 21:38
• Possibly, but it'd need a lot more explanation on SO. – walkytalky Jul 20 '10 at 1:29
• Most programmers know "median". (sort(array))[length/2] is a big enough hint for those who forgot. Also at its most basic for each new point you only need to do a bisection/insert on one half of the array... – Paul Jul 20 '10 at 9:04
• Re-opened following discussion at meta.stats.stackexchange.com/questions/276/… – Rob Hyndman Aug 13 '10 at 0:31
• Too trivial to be more than a comment, but code for the median of 3s is just a + b + c - max(a, b, c) - min(a, b. c). That works fine even if ties are present. That was only obvious to me once I'd thought about it from someone else's code (why is he (in this case) adding and subtracting to get a median???) and a few others might have the same reaction. max() and min() are often implemented as super-fast functions. Sadly there is no such trick in general. – Nick Cox Feb 27 '18 at 15:31

It's bad form to sort an array to compute a median. Medians (and other quantiles) are typically computed using the quickselect algorithm, with $O(n)$ complexity.

You may also want to look at my answer to a recent related question here.

If you're willing to tolerate an approximation, there are other methods. For example, one approximation is a value whose rank is within some (user specified) distance from the true median. For example, the median has (normalized) rank 0.5, and if you specify an error term of 10%, you'd want an answer that has rank between 0.45 and 0.55.

If such an answer is appropriate, then there are many solutions that can work on sliding windows of data. The basic idea is to maintain a sample of the data of a certain size (roughly 1/error term) and compute the median on this sample. It can be shown that with high probability, regardless of the nature of the input, the resulting median satisfies the properties I mentioned above.

Thus, the main question is how to maintain a running sample of the data of a certain size, and there are many approaches for that, including the technique known as reservoir sampling. For example, this paper: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.24.7136

Here is an article describing one possible algorithm. Source code included and a quite serious application (gravitational wave detection based on laser interferometry), so you can expect it to be well tested.

If you maintain a length-k window of data as a sorted doubly linked list then, by means of a binary search (to insert each new element as it gets shifted into the window) and a circular array of pointers (to immediately locate elements that need to be deleted), each shift of the window requires O(log(k)) effort for inserting one element, only O(1) effort for deleting the element shifted out of the window, and only O(1) effort to find the median (because every time one element is inserted or deleted into the list you can update a pointer to the median in O(1) time). The total effort for processing an array of length N therefore is O((n-k)log(k)) <= O(n log(k)). This is better than any of the other methods proposed so far and it is not an approximation, it is exact.

• Could you elaborate on how you propose to do a binary search in a sorted doubly-linked list? – NPE Nov 28 '10 at 12:14
• one 'link' allows you to traverse the list in sorted order; the other allows you to traverse in the order in which the elements appear. It is not clear how you would do this with pointers, though, as @aix questions. – shabbychef Nov 28 '10 at 17:02
• @aix I think your intimation is correct; I would need an indexable skip list, not just a sorted doubly-linked list. The idea is to have a data structure that permits the insertion of one element, the deletion of one element, and finding the median in expected O(log(n)) time (or better). – whuber Nov 28 '10 at 19:17

As you mentioned sorting would be O(n·log n) for a window of length n. Doing this moving adds another l=vectorlength making the total cost O(l·n·log n).

The simplest way to push this is by keeping an ordered list of the last n elements in memory when moving from one window to the next one. As removing/inserting one element from/into an ordered list are both O(n) this would result in costs of O(l·n).

Pseudocode:

l = length(input)
aidvector = sort(input(1:n))
output(i) = aid(n/2)
for i = n+1:l
remove input(i-n) from aidvector
sort aid(n) into aidvector
output(i) = aid(n/2)


Here is a solution O(1) for finding current median, and O(log n) for adding a new number http://www.dsalgo.com/RunningMedian.php

If you can live with an estimate instead of the true median, the Remedian Algorithm (PDF) is one-pass with low storage requirements and well defined accuracy.

The remedian with base b proceeds by computing medians of groups of b observations, and then medians of these medians, until only a single estimate remains. This method merely needs k arrays of size b (where n = b^k)...

I used this RunningStats C++ Library in an embedded application. It is the most simple running stats library I have found yet.