Recalculating or estimating percentile statistics for new data points

Question

After reading The Math Form: Re-Calculating the Standard Deviation I'm curious what sort of statistics can be calculated or estimated in a similar manner.

The post details how—given just the mean, size and standard deviation of a data set—one can recalculate the standard deviation for the same data set plus some new data point without the original data set.

It seems intuitive that there's no possible way to recalculate the median, for example (please correct me if I'm wrong!).

But is there some way to estimate percentiles without retaining the original data set? I'm mostly interested to see what can be done using as little storage as possible, so solutions that rely on indexing are also interesting!

Answer 1

I came across some code that does just what I'm asking when browsing the source for ptaoussanis/tufte, a Clojure performance profiling tool.

ptaoussanis implements a function to merge two sets of statistics, including size, mean, median absolute deviation, and approximations of the 50th, 90th, 95th, and 99th percentiles.

He references Algorithms for calculating variance and a SO question/answer about calculating absolute deviation.

Here's the code snippet, for anyone interested:

(defn merge-stats
  "`(merge-stats (stats c0) (stats c1))` is a basic approximation of `(stats (into c0 c1)))`."
  [m0 m1]
  (if m0
    (if m1
      (let [_ (assert (get m0 :n))
            _ (assert (get m1 :n))

            {^long   n0       :n
             ^long   min0     :min
             ^long   max0     :max
             ^long   sum0     :sum
             ^double mad-sum0 :mad-sum
             ^long   p50-0    :p50
             ^long   p90-0    :p90
             ^long   p95-0    :p95
             ^long   p99-0    :p99} m0

            {^long   n1       :n
             ^long   min1     :min
             ^long   max1     :max
             ^long   sum1     :sum
             ^double mad-sum1 :mad-sum
             ^long   p50-1    :p50
             ^long   p90-1    :p90
             ^long   p95-1    :p95
             ^long   p99-1    :p99} m1

            _ (assert (pos? n0))
            _ (assert (pos? n1))

            n2       (+ n1 n0)
            n0-ratio (/ (double n0) (double n2))
            n1-ratio (/ (double n1) (double n2))

            sum2  (+ sum0 sum1)
            mean2 (/ (double sum2) (double n2))
            min2  (if (< min0 min1) min0 min1)
            max2  (if (> max0 max1) max0 max1)

            ;; Batched "online" MAD calculation here is better= the standard
            ;; Knuth/Welford method, Ref. http://g oo.gl/QLSfOc,
            ;;                            http://g oo.gl/mx5eSK.
            ;;
            ;; Note that there's empirically no advantage in using `mean2` here
            ;; asap, i.e. to reducing (- v1_i mean2).
            mad-sum2 (+ mad-sum0 ^double mad-sum1)

            ;;; These are pretty rough approximations. More sophisticated
            ;;; approaches not worth the extra cost/effort in our case.
            p50-2 (Math/round (+ (* n0-ratio (double p50-0)) (* n1-ratio (double p50-1))))
            p90-2 (Math/round (+ (* n0-ratio (double p90-0)) (* n1-ratio (double p90-1))))
            p95-2 (Math/round (+ (* n0-ratio (double p95-0)) (* n1-ratio (double p95-1))))
            p99-2 (Math/round (+ (* n0-ratio (double p99-0)) (* n1-ratio (double p99-1))))

            mad2 (/ (double mad-sum2) (double n2))]

        {:n n2 :min min2 :max max2 :sum sum2 :mean mean2
         :mad-sum mad-sum2 :mad mad2
         :p50 p50-2 :p90 p90-2 :p95 p95-2 :p99 p99-2})
      m0)
    m1))