Say I have two random streams of two dimensional data. I want to measure how closely their underlying PDF's match. My current method is to estimate the PDF's by accumulating the samples online in a 256x256 histogram. After this, I measure the deviation between the two. I am using sum of absolute differences here, but I believe other metrics will work fine as well.

Some other possibly relevant details: I will be ranking the streams based on their goodness of fit; I don't actually need to perform any hypothesis testing— just find which has highest rank. Also, the distributions are complex, and I don't believe they can be estimated to sufficient accuracy by any well-known analytic distributions.

The problem is, the memory accesses on $t*2^{16}$ bins with $t$ threads, is a significant bottleneck. I was wondering if there were a way to compute the goodness of fit metric without having to accumulate $O(k)$ data about the distributions. Maybe something similar to how there are online calculations of mean/variance that just accumulate one or two variables along the way.

I have looked at Komogorov-Smirnov, Cramer-von Mises, and Shapiro-Wilk tests. They appear to require either 1) an estimation of the sample's CDF, which would still need $O(k)$ variables or 2) ordered data, which isn't possible with an unordered, random stream of data.

Maybe there's some fundamental limitation that prevents a computation comparing PDF's with just $O(1)$ variables?

  • $\begingroup$ 1. If your pmf is really discrete why do you need to accumulate a histogram on the first pass? Its not clear what you're doing there. The KS statistic is based on a difference in cdf, so the question is whether there's a fast update on that. If the pmf is only on a few atoms of probability it may be pretty fast to update. 2. Given a discrete pmf, the C-vM won't require sorting the observations; it should be possible to calculate it from the same cdf you'd update in a K-S, so again the question is whether the cdf is fast to update. $\endgroup$
    – Glen_b
    May 7, 2019 at 23:36
  • $\begingroup$ Maybe I'm misunderstanding the test, but seems you need to first compute the empirical distribution function of the samples. For discrete data, you can calculate that efficiently from a histogram. $\endgroup$
    – Azmisov
    May 7, 2019 at 23:47
  • $\begingroup$ 3. If you're just computing a measure of fit with the S-W (not testing), you don't need to assume normality -- if you have some location-scale family you replace the calculation of expected order statistics. However, I am curious how you're doing it online (since a typical implementation would require something equivalent to sorting the data; however this could be gotten around with a discrete distribution). A chi-square statistic should be very fast to update but would be less sensitive than some other options. $\endgroup$
    – Glen_b
    May 7, 2019 at 23:48
  • $\begingroup$ You just said in your post that a histogram was too inefficient; it can't be both efficient and inefficient. If it is actually efficient, do it. If it is not, stop insisting on doing it that way and instead try to consider if there's a faster way to do it online. It looks to me like the crucial issue is how many values your discrete r.v. can take. If it's say 5 or 6 different values, no problem. If it's thousands, it will be slow to update an unnormalized cdf online. $\endgroup$
    – Glen_b
    May 7, 2019 at 23:50
  • $\begingroup$ Note that if you intend testing any of this you have multiple issues to contend with (such as the impact of the discrete distribution on statistics designed for continuous variables, the sequential testing) $\endgroup$
    – Glen_b
    May 7, 2019 at 23:55

1 Answer 1


Consider updating a chi-squared statistic.

$X^2 = \sum_{i=1}^k \frac{(O_i-np_i)^2}{np_i}$

First let us simplify:

$= \sum_{i=1}^k \frac{O_i^2-2np_iO_i+(np_i)^2}{np_i}$

$= \sum_{i=1}^k \frac{O_i^2}{np_i}-2\sum_{i=1}^k O_i+n\sum_{i=1}^k p_i$

$= \sum_{i=1}^k \frac{O_i^2}{np_i}\,-n$

$= [\frac{1}{n}\sum_{i=1}^k \frac{O_i^2}{p_i}]\,-n$

$= S/n - n$

If we can update $S=\sum_{i=1}^k \frac{O_i^2}{p_i}$ quickly it looks like there's a fast algorithm.

So now let us focus on updating $ S$

Note that $O_i$ is just the count in the $i$th bin. When you add a new observation to bin $i$, you add $(2O_i+1)/p_i$ to the sum - all the contributions from other bins are unchanged, and then you can recompute $S/n-n$ at the new value of $n$, this is $O(1)$ in both $n$ and $k$.


This may potentially be an issue; as $n$ grows you're subtracting two relatively equal quantities (the expected value of the first term $S/n$ should be about $n+k$(*) if the null is true but will tend to be larger otherwise). As $n$ grows very large relative to $k$ this may become an issue; specifically, computing something akin to $(1+k/n) - 1$ may be inaccurate if $k/n$ is very small.

At the same time if $n$ is really large, $S$ may be large relative to the increment in $S$. Ideally we avoid both issues by adjusting the S-calculation.

I believe the same sort of idea used in online variance update calculations could be adapted to this computation in that case, it would slightly complicate the update of $S$ - in a way that it's already reduced by $n^2$ (or something close to it).

Actually I think we just add something like $(2O_i+1)/p_i - 2n-1$ to $S$ instead and don't subtract $n$ at the final step (or it may be $-2n+1$ depending on when we update $n$).

(However if $n$ is only a few million it's probably not worth it; if you have several figures of accuracy on your statistic - and you should have plenty to spare - it will probably be sufficient to your purposes. If $n$ is many billions, you should do it but you will want to consider efficiency of the update pretty carefully then)

* actually, n+k-1 but that matters not in the least to the issue.

  • $\begingroup$ So maybe this is my fault for not clarifying precisely, and I apologize. I meant an online algorithm in that it reads the data one-by-one in a single pass. I don't actually need to compute the final metric on each iteration. In my original attempt, I'd accumulate $O(k)$ bins in a single pass, then compute the metric. The binning being the bottleneck, I was wondering if you could do it with just $O(1)$ memory total, over all iterations. $\endgroup$
    – Azmisov
    May 8, 2019 at 2:31
  • $\begingroup$ My line of thinking was something like: Assume $p_i$ is already known (accumulated before for one of the data streams). For each observation you see increment your metric by $p_i/n$, where $p_i$ is the observation's corresponding bin. That gives you $\sum p_io_i$, where $o_i$ is probability from 2nd data stream's bin... so maybe you could compare this to $\sum p_i^2$? $\endgroup$
    – Azmisov
    May 8, 2019 at 2:31
  • $\begingroup$ If your $p_i$ are themselves estimates of true p's based on data, the entire question is incorrectly specified and you should write a new one that explains that you want to compare proportions across data streams (or whatever it is you want). But if you have $p$s prespecified (as the question here suggests) then I don't see what you're getting at with this at all. In any case I don't see what you mean by "increment your metric" (what metric?) nor "where $p_i$ is the observations corresponding bin" (no, $p$ is a probability, $i$ is the bin, $n$ is the index of the current observation). $\endgroup$
    – Glen_b
    May 8, 2019 at 3:05
  • $\begingroup$ If you're trying to do it without storing either the $p_i$s or the $O_i$s that's going to be tricky. If you do store $p$, storing $O$ doesn't change the big-$\mathcal{O}$ factor on storage space, just the multiplier out the front. Anyway, I am presently confused (and I presumably missed something important), but maybe you can clarify. $\endgroup$
    – Glen_b
    May 8, 2019 at 4:04
  • $\begingroup$ I was coming from the angle, "if its impossible to compute in $O(1)$ when one of the distributions is known (e.g. $p_i$ known apriori), then its probably impossible when both are unknown." So it was just an intuition pump that got me thinking it might actually be possible. To clarify you've got $o_i$=1st distribution's histogram probabilities, $p_i$=2nd distribution's histogram probabilities (assume you estimated it already— actually unknown in real problem). Then if you make $N$ observations, you can compute $S=\sum_{n}^{N}p_i/N$, where $i$ is sample $n$'s bin. Giving you $S=\sum_i p_io_i$. $\endgroup$
    – Azmisov
    May 8, 2019 at 4:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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