Efficient computation of multivariate empirical distribution functions at the observed values

Consider the evaluation of model-based functions of cumulative distribution functions that are integrals. When the cumulative distribution function does not have a tractable form but simulation of the multivariate distribution is easily feasible, we can evaluate the integral via a Monte Carlo sample, replacing the model-based distribution function by the empirical distribution function. Given a simulation sample of size N, the naive method uses $$O(N^{2})$$ O ( N 2 ) comparisons to compute the empirical distribution function at all N sample vectors. To obtain faster computational speed when N needs to be large to achieve a desired accuracy, we propose methods modified from the popular merge sort and quicksort algorithms that preserve their average $$O(N\log _{2}N)$$ O ( N log 2 N ) complexity in the bivariate case. The modified merge sort algorithm can be extended to the computation of a d-dimensional empirical distribution function at the observed values with $$O(N\log _{2}^{d-1}N)$$ O ( N log 2 d - 1 N ) complexity. Simulation studies suggest that the proposed algorithms provide substantial time savings when N is large.