I am trying to calculate the "overall" median of a variable that is spread across two datasets. I have access to the raw data in each dataset but can't bring their raw data together. What would be the best way to calculate the median (or any quantile for that matter)? What I was thinking to do was to get the quantiles from the "remote" dataset, say with N_1 samples and bring them to the other dataset; then sample N_1 values from these quantiles (essentially discrete empirical CDF) and use these values as "raw" scores to calculate the overall median or any other quantile. I wonder if there is a more sophisticated way of doing this?

I can't combine the two datasets because of data governance issues, where I have to analyse them separately on different computers.

  • 2
    $\begingroup$ Could you explain what you mean by "can't bring the raw data together"? $\endgroup$ – whuber Nov 11 '20 at 20:04
  • 1
    $\begingroup$ I've edited the questions. Thanks. $\endgroup$ – A. Blizzard Nov 11 '20 at 22:35
  • 1
    $\begingroup$ Thank you. See this search. stats.stackexchange.com/questions/7959 is directly applicable. It would help to know (a) how accurately you need to compute the median, (b) what information you can carry over from one system to the other, and (c) how often you can iterate back and forth between systems. $\endgroup$ – whuber Nov 11 '20 at 22:39
  • 1
    $\begingroup$ By means of a range of quantiles (along with the dataset sizes) you can closely approximate the distribution in each system; the combined dataset is the mixture of the two, whose quantiles are easily calculated. As an example of how you could go about choosing suitable quantiles, see stats.stackexchange.com/a/35268/919 for one proposal that provides good control over the precision. $\endgroup$ – whuber Nov 11 '20 at 23:10
  • 1
    $\begingroup$ OK, thanks! I think mixture is the keyword here. $\endgroup$ – A. Blizzard Nov 11 '20 at 23:20

By adapting the (simple, recursive) method for approximating QQ plots I described at https://stats.stackexchange.com/a/35268/919, you can efficiently find a precise (piecewise linear) approximation of the full empirical quantile function using just a few values, typically a few hundred to a few thousand even for a batch of a billion numbers (depending on how precise the summary needs to be). It is fast and straightforward to combine two such summaries to produce the full empirical quantile function of the combined datasets.

This algorithm requires $O(n\log(n))$ computational effort to process a batch of $n$ numbers and, because it is a divide-and-conquer algorithm, it is highly parallelizable. This makes it practicable even for large datasets; for instance, on a single core R requires only about four minutes to process a dataset of one billion ($10^9$) numbers.

The crux of the matter is to compute the quantile function $F^{-1}$ of a mixture of two distributions $F_X$ and $F_Y$ in terms of those distributions and their quantile functions $F_X^{-1}$ and $F_Y^{-1}.$ Let $p_X$ and $p_Y,$ with $p_X+p_Y=1,$ be the (positive) mixture weights. This means that for any number $x,$

$$F(x) = p_X F_X(x) + p_Y F_Y(x).\tag{*}$$

Let $0 \lt q \lt 1$ (the question concerns the median, where $q=1/2$) and suppose we want to find the quantile $F^{-1}(q).$ By definition, this is the smallest $x$ among the numbers for which $F(x) \ge q.$ Replacing $F$ by the foregoing expression and algebraically rearranging gives

$$F_X(x) \ge \frac{1}{p_X}\left(q - p_Y F_Y(x)\right),$$

which means

$$x = F_X^{-1}\left(\frac{1}{p_X}\left(q - p_Y F_Y(x)\right)\right).$$

A simple search for the smallest zero of

$$g(x) = F_X^{-1}\left(\frac{1}{p_X}\left(q - p_Y F_Y(x)\right)\right)-x\tag{**}$$

thereby computes the $q$ quantile of $F.$

As an example, here are two small datasets of 700 and 300 numbers each, named x and y.

Figure 1: histograms of x, y, and the combined data.

The challenge is to approximate the median (or any other quantile) of the left (combined) distribution based on small summaries of its component distributions (shown in color). For these data, the proposed algorithm summarizes x using $34$ (value, quantile) pairs and it summarizes y using just $10$ (value, quantile) pairs. These summaries approximate the quantile functions to within a maximum error of $0.5\%$ of the range of all the data. (The amount of allowable error can be set to anything, even zero. Zero is an effective choice for discrete data: in this case, the algorithm is equivalent to tabulating the values in each dataset.)

Here is a plot of the empirical distribution function of all the data, $F$ (black), on which is superimposed the mixture distribution $(*)$ as computed from the approximate ECDFs for $F_X$ and $F_Y$ (shown in red):

Figure 2

There are two vertical lines shown: a gray one is located at the true median of the combined data while the red one is located at the root of $g$ (equation $(**)$) computed using the approximate ECDFs. Indeed, they differ by only $0.03,$ which is less than $0.5\%$ of the full range of the data, as intended.

This is the R code that generated the examples. It implements a quantile/distribution function approximator, quantilefun, and shows how to use it for computing quantiles of mixtures of two distributions. It uses the built-in, robust uniroot function to find the mixture quantiles. Because the approximating functions are piecewise linear, they are fast to evaluate, making the search for roots reasonably efficient.

# Compute an approximation to the empirical quantile function.
# It will be accurate everywhere to within the given relative tolerance
# (a multiple of the range of `x`).
quantilefun <- function(x, tol=0.005) {
  x <- sort(x[!is.na(x)])
  n <- length(x)
  if (n <= 1) stop("need at least two non-NA values to interpolate.")
  threshold <- tol * diff(range(x))
  # Find an index i <= k <= j at which x[k] departs most from the linear
  # interpolator of x[i] and x[j].  Assumes i <= j.
  # Returns a list of the index and the discrepancy.
  f <- function(i,j) {
    beta <- (x[j] - x[i]) / (j - i) 
    dy <- abs(x[i:j] - x[i] - beta * (i:j - i))
    k <- which.max(abs(dy))
    list(k=k+i-1, dy=dy[k])
  # Locate the key indexes into `x` between `i` and `j`.
  # Returns them in ascending order.
  keys <- function(i, j) {
    k <- f(i, j)
    if (abs(k$dy) <= threshold) return(integer())
    c(keys(i,k$k), k$k, keys(k$k,j))
  k <- c(1, keys(1, n), n)
  q <- (k-1)/(n-1)
  # Although only `Indexes` and `Values` are needed, additional summaries
  # are returned for ease of use.
       Quantile=approxfun(q, x[k], yleft=x[1], yright=x[n]), 
       F=approxfun(x[k], q, yleft=0, yright=1))
# Create test data.
# It is most difficult to combine datasets with very different quantiles.
x <- rgamma(7e2, 2, 2/3)
y <- rgamma(3e2, 4, 2) + 5
# Compute the approximate empirical distribution functions.
tol <- 5e-3
r.x <- range(x)
r.y <- range(y)
n.x <- length(x)
n.y <- length(y)
r <- range(c(r.x, r.y))
n <- n.x + n.y
  f.x <- quantilefun(x, tol * diff(r) / diff(r.x))
  f.y <- quantilefun(y, tol * diff(r) / diff(r.y))
# Estimate the median of the combined data using only the approximate ECDFs.
  f <- function(x, q=1/2) f.x$Quantile(n/n.x * (q - n.y/n * f.y$F(x))) - x
  median.combined <- uniroot(f, r, tol=1e-12)$root
# Display the median and the estimated median.
median <- median(c(x, y))
d <- -floor(log10(f.x$Threshold))
   `Combined estimate`=round(median.combined, d),
   Threshold=round(f.x$Threshold, d)))
(c(N = length(f.x$Indexes) + length(f.y$Indexes)))
# Plot the combined ECDF.
if (n < 1e5) {
  plot(ecdf(c(x,y)), main="ECDF of all Data", ylab=expression(F(x)))
  # Superimpose the combined approximations.
  curve((n.x * f.x$F(x) + n.y * f.y$F(x)) / n, add=TRUE, col="#f0000080", lwd=2)
  # Plot the median and its estimate.
  abline(v = c(median, median.combined), col=c("Gray", "Red"), lwd=2)
  abline(h = 1/2, lwd=2)
# Plot the histograms.
hist(c(x,y), freq=FALSE, col="Gray")
hist(x, xlim=range(c(x,y)), freq=FALSE, col="#ffd0d0")
hist(y, xlim=range(c(x,y)), freq=FALSE, col="#d0d0ff")
  • $\begingroup$ whuber thank you very much for your help and your time. Never new about uniroot function in R, new only about optim, I was actually doing numerical search on the mixture function *. I will adapt this to arbitrary number of datasets. I can't export exact values, I'll have to jitter them a bit. $\endgroup$ – A. Blizzard Nov 12 '20 at 18:01
  • $\begingroup$ Jittering is a good idea. If you jitter with a maximum error of $\epsilon/2$ and then approximate the ECDF of the jittered data with a maximum error of $\epsilon/2,$ you will still control the overall error of the result to within $\epsilon.$ Dealing with mixtures of more than two components starts getting computationally expensive. You might have to mix two components and then approximate the approximation of their mixture by computing a set of predetermined quantiles (such as all percentiles) using uniroot, then iterate that process over the remaining components. $\endgroup$ – whuber Nov 12 '20 at 18:06
  • $\begingroup$ OK. Thank you so much again! $\endgroup$ – A. Blizzard Nov 12 '20 at 18:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.