# Algorithms for computing multivariate Empirical distribution function (ECDF)?

One dimensional ECDF is fairly easy to compute. When it comes to two dimensions and up, however, online resources become sparse and hard to reach. Can anyone suggest, define and/or present efficient algorithms (not ready made implementation) for computing multivariate ECDF?

• This may well be a computer science question, but I think this is the best place to find an answer, let me know if I should look elsewhere. Thanks. Commented Aug 2, 2016 at 18:40
• Is there really any fundamental difference? Computing a univariate ECDF is equivalent to sorting the data. Computing a multivariate ECDF is equivalent to sorting the data lexicographically.
– whuber
Commented Aug 2, 2016 at 19:28
• @whuber, not exactly, as far as I know. For each data-point, X(i), we need to count the number of points contained in the hypercube that is defined by it (from -inf up to and including X(i) in all dimensions). Lexicographical (dictionary?) sorting will not necessarily work here, as the data-points have to be compared in every dimension separately. E.g.: (2,3,4) will be lexicographically greater compared with (1,2,15), but the hypercube that's defined by (2,3,4) will not contain (1,2,15) since 15>4. Commented Aug 2, 2016 at 20:54
• True, the correspondence is not that direct. But one would exploit that sorting, or something like it, to build a point quadtree (or octree etc) with $O(n\log(n))$ effort. You might want to investigate the computational geometry and spatial indexing literature for details.
– whuber
Commented Aug 2, 2016 at 22:22

On further investigation, the following paper gives efficient algorithms for the k-D ECDF problem:

Bentley, J. L. (1980). Multidimensional divide-and-conquer. Communications of the ACM, 23(4), 214-229.

The main data structure introduced is known as a range tree, and is somewhat similar to a k-d tree, but uses a space-for-time tradeoff to achieve faster range queries. The author of the above paper, Jon Bentley (of Programming Pearls fame), is the inventor of both data structures.

Both are binary trees which recursively partition a set of $k$ dimensional points by splitting along a coordinate axis at the median. In a k-d tree the sub-trees of a node are split along the $d$-th dimension, where $d$ cycles through $1\ldots k$ moving down the tree. In a range tree the sub-trees are always split along the first dimension, but each node is augmented with a $k-1$ dimensional range tree defined over the remaining dimensions.

At the time of this writing, the Wikipedia page for "Range Tree" linked above points to a CS lecture (Utrecht U.) comparing these two tree types from circa 2012. This suggests that these data structures are still essentially "state of the art". There is mention of an improved "fractional cascading" variant for range trees, but for the all-points ECDF problem this just allows the performance of Bentley's algorithm to be achieved via repeated queries of the range tree.

• Thanks for the interesting paper! I think this is what I need re: k-d trees. Would be great to see alternative methods. Unless this one is the state of the art. Commented Aug 25, 2016 at 15:50
• @AlexanderF. I updated my answer to better describe the algorithm (including a more "official" reference). It appears the approach is still close to state of the art. For recent developments the key phrase appears to be "orthogonal range queries" if you want to investigate further. Commented Aug 27, 2016 at 2:57

I am not sure if there is a more efficient way to compute the ECDF at the data points, but the following brute force approach should be efficient for computing the ECDF over the data "grid". It is a simple generalization of the 1D version.

Assume you have a data set consisting of $N$ points in $d$ dimensions, given in the $N\times d$ matrix $X$. For simplicity I will assume that $X$ consists entirely of unique numbers (i.e. general position*). I will use Matlab notation in the following pseudo-code, as it is how I thought of the algorithm, but I can expand on this if there is interest.

First compute

$[x_{:,k},I_{:,k}]=\text{sort}[X_{:,k}]$ for $k=1:d$,

where $I$ is the coordinate-wise rank matrix, and $x$ is the coordinate-grid axis matrix (both of size $N\times d$).

Then rasterize the data points into the implied data grid, computing an (normalized) histogram as $P=\text{accumarray}[I,\frac{1}{N},N\times\text{ones[1,d]}]$.

Then integrate this "EPDF" in each dimension to give the ECDF: $P=\text{cumsum}[P,k]$ for $k=1:d$.

Now $P_{i_1,\ldots,i_d}$ is the ECDF sampled at $x_{i_1,1},\ldots x_{i_d,d}$.

This algorithm takes time $\text{O}[N\log N]$ for each sort and $\text{O}[N^d]$ for each sum, so the total cost is $\text{O}[d(N^d+N\log N)]$. As the gridded ECDF itself has $\text{O}[N^d]$ elements, this should be essentially optimal.

(*The assumption of distinct points can be relaxed by using $\text{unique}[]$ instead of $\text{sort}[]$, along with a bit of book-keeping.)

• You might be interested in learning about quadtrees and their higher-dimensional generalizations, which provide efficient ways to search Euclidean spaces for points. They asymptotically use $O(N\log(N))$ resources, which is far better than $O(N^d)$ for $d\gt 1$.
– whuber
Commented Aug 24, 2016 at 15:47
• @whuber I have some idea of these, (e.g. k-d trees). I am not sure if there is a single "best answer" here? Typically for a problem like this, you will also specify what operations your abstract ECDF data structure should support (e.g. point queries, subspace integrals, updating with new points, etc.). This will help determine what implementation is best suited. Commented Aug 24, 2016 at 15:53
• I believe it should be clear what operations need to be supported for an ECDF. The minimal one is to evaluate it at any point in the space. It is true that if one is going to build an ECDF dynamically, alternative approaches might be superior, but those issues appear to be beyond the scope of the present question.
– whuber
Commented Aug 24, 2016 at 15:56
• @GeoMatt22, this indeed looks like a method for computing histogram and it may be fine in cases where approximation is "good enough". However, why use a method which is O(N^d) when the brute force approach in O(d*N^2). Fro example, for now I don't have too large of a data-set so I use the following Matlab one liner to compute d-dimensional ECDF with O(d*N) storage complexity (C(i) is the frequency of data-point Y(i,:)): arrayfun(@(i) sum(C(all(bsxfun(@le,Y, Y(i,:)), 2))), (1:size(Y,1)).'); Commented Aug 25, 2016 at 15:47
• (+1) Not for giving an efficient algorithm, but for clearly explaining an inefficient one that helped me understand the problem. Commented Sep 7, 2016 at 8:13