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?
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.
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.
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.
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
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 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 size() const {return tot;};
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;
}
}
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.
Here's an answer to the question asked on stackoverflow: https://stackoverflow.com/questions/1058813/on-line-iterator-algorithms-for-estimating-statistical-median-mode-skewness/2144754#2144754
The iterative update median += eta * sgn(sample - median) sounds like it could be a way to go.
