Algorithms for computing multivariate Empirical distribution function (ECDF)?

Question

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?

Answer 1

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.

Answer 2

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