# What is a good algorithm for estimating the median of a huge read-once data set? [duplicate]

I'm looking for a good algorithm (meaning minimal computation, minimal storage requirements) to estimate the median of a data set that is too large to store, such that each value can only be read once (unless you explicitly store that value). There are no bounds on the data that can be assumed.

Approximations are fine, as long as the accuracy is known.

Any pointers?

• You can try to find a median based on grouped frequency distribution, here is some details
– ryfm
Jul 20, 2010 at 20:12
– user28
Jul 21, 2010 at 15:33
• @Srikant:> it's a pretty active area of research in statistics :) The solution closest to the lower theoretical bounds in terms of storage involve some pretty clever probability constructs as well. All in all i was surprised when i first looked unto it a couple of month ago; there is more stats here than meets the eye. Oct 9, 2010 at 14:21

Could you group the data set into much smaller data sets (say 100 or 1000 or 10,000 data points) If you then calculated the median of each of the groups. If you did this with enough data sets you could plot something like the average of the results of each of the smaller sets and this woul, by running enough smaller data sets converge to an 'average' solution.

• This is interesting, and where some statistical advice could come in! Assume in total I've got (say) 500,000 i.i.d. points and I look at groups of (say) 1,000 of them, and calculate the median of each group. Now I've got 500 medians. Is there theory that could allow me to calculate a confidence interval for the overall median based on these 500 medians? Jul 21, 2010 at 17:10
• So, according to a long lost colleague, the best apropoach seems to be Chiranjeeb Buragohain and Subhash Suri. Quantiles on Streams. cs.ucsb.edu/~suri/psdir/ency.pdf I also like Ian's approach, as these medians of smaller data sets will converge to a normal distribution, and so I can form conf intervals for the medians. Jul 27, 2010 at 15:10

How about something like a binning procedure? Assume (for illustration purposes) that you know that the values are between 1 and 1 million. Set up N bins, of size S. So if S=10000, you'd have 100 bins, corresponding to values [1:10000, 10001:20000, ... , 990001:1000000]

Then, step through the values. Instead of storing each value, just increment the counter in the appropriate bin. Using the midpoint of each bin as an estimate, you can make a reasonable approximation of the median. You can scale this to as fine or coarse of a resolution as you want by changing the size of the bins. You're limited only by how much memory you have.

Since you don't know how big your values may get, just pick a bin size large enough that you aren't likely to run out of memory, using some quick back-of-the-envelope calculations. You might also store the bins sparsely, such that you only add a bin if it contains a value.

Edit:

The link ryfm provides gives an example of doing this, with the additional step of using the cumulative percentages to more accurately estimate the point within the median bin, instead of just using midpoints. This is a nice improvement.

• The problem with the binning approach is that we do not have a good upper bound for the data, and so the midpoint for the largest bin would have to be huge. So, we'd need a huge number of bins (not enough memory for that), or have pretty wide bins (which would then lead to a fairly inaccurate answer.) And the data is not very sparse. Jul 21, 2010 at 15:13
• Since you are only interested in the median why couldn't you make the bins wider at higher values of your variable? Jul 27, 2010 at 0:23
• drknexus - because we don't know what the largest bin should be. Jul 27, 2010 at 14:58
• Do you have any intuition as to what the range will be? If you're fairly sure that over half of the answers will be below number N, then you can make your last bin as large as you want. Maybe your last bin is all numbers greater than 1 trillion - would that be high enough? With the amount of memory in modern systems you can store a LOT of bins and achieve fairly high resolution. In terms of data structures, we're not talking anything fancy and memory intensive here. Jul 27, 2010 at 22:17
• Any intuition? yes. And your approach could work in general. However, in this case we can not have a lot of memory/computation. It is in a networking application where the device could see tens of thousands of items per second, and have VERY little processing left over for this purpose. Not the ideal/typical scenario, I know, but that is what makes it interesting! Jul 28, 2010 at 16:43

I re-direct you to my answer to a similar question. In a nutshell, it's a read once, 'on the fly' algorithm with $O(n)$ worst case complexity to compute the (exact) median.

I implemented the P-Square Algorithm for Dynamic Calculation of Quantiles and Histograms without Storing Observations in a neat Python module I wrote called LiveStats. It should solve your problem quite effectively.

• (+1) Thanks for stopping by and providing that link, Sean!
– whuber
Mar 9, 2013 at 22:23

The Rivest-Tarjan-Selection algorithm (sometimes also called the median-of-medians algorithm) will let you compute the median element in linear-time without any sorting. For large data sets this is can be quite a bit faster than log-linear sorting. However, it won't solve your memory storage problem.

I've never had to do this, so this is just a suggestion.

I see two (other) possibilities.

Half data

1. Load in half the data and sort
2. Next read in the remaining values and compare against the your sorted list.
1. If the new value is larger, discard it.
2. else put the value in the sorted list and removing the largest value from that list.

Sampling distribution

The other option, is to use an approximation involving the sampling distribution. If your data is Normal, then the standard error for moderate n is:

1.253 * sd / sqrt(n)

To determine the size of n that you would be happy with, I ran a quick Monte-Carlo simulation in R

n = 10000
outside.ci.uni = 0
outside.ci.nor = 0
N=1000
for(i in 1:N){
#Theoretical median is 0
uni = runif(n, -10, 10)
nor  = rnorm(n, 0, 10)

if(abs(median(uni)) > 1.96*1.253*sd(uni)/sqrt(n))
outside.ci.uni = outside.ci.uni + 1

if(abs(median(nor)) > 1.96*1.253*sd(nor)/sqrt(n))
outside.ci.nor = outside.ci.nor + 1
}

outside.ci.uni/N
outside.ci.nor/N


For n=10000, 15% of the uniform median estimates were outside the CI.

• The data set is potentially too big to read in half of it...it is in a networking context where the device doing the processing can see tens of thousands of items per second, and probably has enough memory to store only a few hundred. Also the data is definitely not Gaussian. In fact it does not fit well to any of the common distributions. Jul 21, 2010 at 15:17

The Remedian Algorithm (PDF) gives a one-pass median estimate 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)...

If the values you are using are within a certain range, say 1 to 100000, you can efficiently compute the median on an extremely large number of values (say, trillions of entries), with an integer bucket (this code taken from BSD licensed ea-utils/sam-stats.cpp)

class ibucket {
public:
int tot;
vector<int> dat;
ibucket(int max) {dat.resize(max+1);tot=0;}

int operator[] (int n) const {
assert(n < size());
int i;
for (i=0;i<dat.size();++i) {
if (n < dat[i]) {
return i;
}
n-=dat[i];
}
}

void push(int v) {
assert(v<dat.size());
++dat[v];
++tot;
}
};

template <class vtype>
double quantile(const vtype &vec, double p) {
int l = vec.size();
if (!l) return 0;
double t = ((double)l-1)*p;
int it = (int) t;
int v=vec[it];
if (t > (double)it) {
return (v + (t-it) * (vec[it+1] - v));
} else {
return v;
}
}

• Also, this can be extended to using a finite number of bins for real-time medians, etc. Aug 20, 2014 at 15:27

Another thought is in the line of random sampling. I had similar problem. My problem is that I have > 100 million data (each has 10 million). The computation took too long. You can just random sample N data points from the 10 million, then find the median on those.

• How does the efficiency of this procedure compare to the other proposed solutions in the other answers?
– Sycorax
Oct 30, 2020 at 1:52
• It is not about efficiency, it is about lazyness. In my case, I know the median is between 1 and 151. I got an huge array of > 2 million each. I know by sampling just a few thousands, the median will be very close the the median obtained from the original. This is a technique used frequently in statistics. Oct 30, 2020 at 4:12
• Some data structures are based on probability. Oct 30, 2020 at 4:13
• Yes, I've heard of random sampling before. You might be surprised to learn that in statistics, the word "efficiency" characterizes the precision of a result in terms of the sample size. In other words, the efficiency of an estimator provides an exact answer to the question "how close will I be to the median obtained from the original?"
– Sycorax
Oct 30, 2020 at 4:21
• Can you characterize that trade-off more precisely? How much loss in precision are you trading for the increase in speed?
– Sycorax
Oct 30, 2020 at 4:28