Estimating mean and variance / standard deviation

[npiguet]

I was wondering if there was a way to calculate the mean and standard deviation of the distribution summarized by a T-Digest?

Disclaimer: I don't really understand all the clever math behind T-Digest

Mean:
It seems that calculating the average of all centroids weighted by their respective count should get me close enough.

Variance:

I was thinking of:

1. Calculate the mean
2. For each centroid, calculate (centroid-mean)^2 * count, and then sum the result from all centroids
3. Then the variance is equal to the value calculated at point 2, divided by sum(centroid.count)

This is obviously wrong, since it makes the assumption that all the sample represented by a centroid have the same value, but I can't really figure out a reasonable way to get more accurate values.

Do you know a better way to calculate the mean and variance? Would it be possible to add the corresponding methods to the TDigest class?

[tdunning]

Your intuition is right on.

The mean needs to be a weighted mean, of course. But other than that, what you said was fine.

The computation of the variance is fine. We purposely have lost the information about location of samples contributing to a centroid, but we do know that the mean is correct (because we only merge centroids, never split them). Because the variance is summed squares, we can treat them as if they are concentrated at the mean. That's pretty handy.

So, let's write it down somewhat formally. Lets say that sample $i$ is $x_i$. The samples are grouped into sets which contribute to a particular centrod where $s_j$ is the $j$-th set of indices for the samples in that set. That let's us write per set means and combine them:

$$\mu_j = \frac 1 {N_j} \sum_{i \in s_j} x_i$$

$$\mu = \frac 1 {N} \sum_j N_j \mu_j$$

Then we can use a clever trick that comes from the variance being quadratic. The $j$-th variance is

$$\sigma_j^2 = \frac 1 {N_j} \sum_{i \in s_j} (x_i - \mu_j)^2 = \frac 1 {N_j} \left( \sum_{i \in s_j} x_i^2 - 2\mu_j x_i + \mu_j^2\right)$$

$$= \frac 1 {N_j} \left( \sum_{i \in s_j} x_i^2 - 2\mu_j \sum_j x_i + N_j \mu_j^2 \right)$$ $$= \frac 1 {N_j} \sum_{i \in s_j} x_i^2 - \mu_j^2 $$

This is cool since we can use it to relate $\sum (x_i - \mu_j)^2$ to $\sum (x_i - \mu)^2$. Thus we have this:

$$\sigma^2 = \frac 1 {N} \left[ \sum_j \sum_{i\in s_j} (x_i - \mu^2) \right]$$ $$ = \frac 1 {N} \left[ \sum_j N_j \left(\sigma_j^2 + \sum_j (\mu_j^2 - \mu^2)\right) \right]$$

The second part of this is what you had. The first part depends on the individual variances for each centroid which we don't really know (although we could retain those values if we wanted to double the size of the t-digest). We could assume that these are all really small. This is kind of right. But I think we can do better.

For every centroid $j$ where $N_j = 0$, the centroid variance really is zero. Since the first and last centroids are guaranteed to have only a single sample, that handles that. For other centroids that have more than one sample, we can get a better approximation by assuming that the samples for the centroid are uniformly distributed between the boundary points with adjacent centroids. Let the boundary between centroid $j$ and centroid $j+1$ be $b_j$. Then we can have a good approximation of $\sigma_j$ for every non-trivial centroid as

$$\sigma_j^2 \approx (b_{j+1}-b_j)^2/12$$

Where $j \in \lbrace 2 \ldots n-1 \rbrace$ where $n$ is the number of centroids.

Does this make sense? (TLDR, you are right but with a small potential for improvement)

[npiguet]

Thanks for the reply. Since the code for MergingDigest does a linear interpolation while calculating quantiles that fall somewhere in the middle of a bucket, I did have the intuition that the variance of that linear distribution should play a role as well. I wasn't sure what that was, but your post above made it clear.

Now I need to be able to relate this to the code. Can you confirm that we have:

$\mu_j = centroids[j].centroid$
$N_j = centroids[j].count$

Also, I can't find a clear reference to the boundaries between centroids in the code. Can derive the boundary between centroids from the List<Centroid> returned by TDigest.getCentroids() using something like:

$b_j=\frac{N_j \mu_j + N_{j+1} \mu_{j+1} }{N_j + N_{j+1}}$ ?

PS: I'm sure it's correct, but I couldn't quite follow the development in your second set of equations beyond the first line :-)

[tdunning]

That definition of $b_j$ is just what I had in mind except I think you need to reverse $N_j$ and $N_{j+1}$. And your interpretations of $\mu_j$ and $N_j$ are also correct.

Thus

$$b_j = \frac {N_{j+1} \mu_j + N_j \mu_{j+1}} {N_j + N_{j+1}}$$

The idea is that $b_j$ should be further from $\mu_j$ if $N_j$ is larger.