2
$\begingroup$

I'm trying to re-implement the Ckmeans.1D.DP algorithm in Python. I have the actual dynamic optimization part down, but I'm a little confused by the BIC computation they use for selecting the number of clusters. I'm not concerned right now about whether this is the best technique; I'm just trying to achieve parity with their implementation.

In particular, this variance computation (performed inside each cluster) is baffling:

double mean = 0.0;
double variance = 0.0;

for (size_t i = indexLeft; i <= indexRight; ++i) {
  mean += x[i];
  variance += x[i] * x[i];
}
mean /= numPointsInBin;

if (numPointsInBin > 1) {
  variance = (variance - numPointsInBin * mean * mean) / (numPointsInBin - 1);
} else {
  variance = 0;
}

(source)

I read this as follows: $$\begin{align} \hat \mu_k &= \frac{1}{n_k} \sum_{i=i_\text{left}}^{i_\text{right}} x_i \\ \hat \sigma^2_k &= \begin{cases} \frac{1}{n_k - 1} \sum_{i=1_\text{left}}^{i_\text{right}} \left( x_i^2 - \hat \mu_k^2 \right) &,\ n_k > 1 \\ 0 &,\ n_k = 1 \end{cases} \end{align}$$

where $\hat \mu_k$ is a cluster mean, $\hat \sigma^2_k$ is a cluster variance, $n_k$ is a cluster size, and $k$ indexes the clusters.

This is totally bizarre: the maximum likelihood estimator for variance is $\frac{1}{n - 1} \sum_i \left(x_i - \hat \mu\right)^2$ (for generic $n$, $i$, and $\hat \mu$), and it's pretty clear that $\left(x_i - \hat \mu\right)^2 = x_i^2 - 2 x_i \hat\mu + \hat\mu^2 \neq xi^2 - \hat\mu^2$. What gives?

Improve this question
$\endgroup$
10
  • $\begingroup$ This is just the standard garden variety textbook one-pass formula for computing sample variance with a denominator of (n-1); with the meaning of the variance variable changing along the way to make it a little more confusing..Note that MLE estimator actually has n, not n-1, in the denominator. On your "pretty clear" thing in your penultimate sentence, if you take expectations, you'll see it comes out o.k. This is a numerically unstable (potentially inaccurate) way of calculating sample variance, but would be fine if all calculations were done in exact arithmetic, i.e. with no roundoff error. $\endgroup$
    – Mark L. Stone
    Commented May 2, 2016 at 17:08
  • $\begingroup$ Although this result is difficult to search for, it has been shown in dozens of threads. The first relevant hit turned up by that search is stats.stackexchange.com/questions/97746, where a generalization (to regression involving more than just a constant, as here) is proven in several ways. But since you don't recognize the formula for $\hat\sigma^2_k$ as the standard definition of a sample variance, exactly what do you understand the variance to be? $\endgroup$
    – whuber ♦
    Commented May 2, 2016 at 17:18
  • 2
    $\begingroup$ Given that this numerically unstable code is used in an archival submission (CRAN), it is quite disconcerting, and contributes to the evidence that you can't automatically trust anything in CRAN. The package authors should learn basic numerical analysis. They would then realize how unsuitable their code is for production use. They could have used a numerically stable one-pass algorithm. ... Although it least this isn't being done in single precision on a GPU with billions of data points, as some folks do. $\endgroup$
    – Mark L. Stone
    Commented May 2, 2016 at 17:20
  • $\begingroup$ @whuber I figured it out, I just went out to lunch before I posted my answer :) $\endgroup$
    – shadowtalker
    Commented May 2, 2016 at 17:31
  • 2
    $\begingroup$ Yes, that's correct. As far as the claims of stability go, it comes down to something I'm sure you know. If you wanted (say) to compute the variance of the data $1,2,3$, you could do so accurately. You also know it's the same as the variance of the data $1000000001,1000000002, 1000000003$. However, the squares of all those numbers (as double precision floats) are large and subtracting equally large numbers from them can lose all the precision inherent in the format. It happens all the time when clustering geographic coordinates! $\endgroup$
    – whuber ♦
    Commented May 2, 2016 at 19:17

1 Answer 1

Reset to default
1
$\begingroup$

The answer (which I discovered in the process of writing the question) is to look at the forest, not the trees!

Expand the sum: $$\begin{align} \sum_i \left( x_i - \hat\mu \right)^2 &= \sum_i \left( x_i^2 - 2\hat\mu x_i + \hat\mu^2 \right) \\ &= \sum_i x_i^2 - 2\hat\mu\sum_i x_i + n\hat\mu^2 \end{align}$$

Then the trick is to realize that $\sum_i x_i = \frac{n}{n}\sum_i x_i = n\hat\mu$:

$$\begin{align} \sum_i \left( x_i^2 \right) - 2\hat\mu \cdot n\hat\mu + n\hat\mu^2 &= \sum_i \left( x_i^2 \right) - 2\hat\mu^2 + n\hat\mu^2 \\ &= \sum_i \left( x_i^2 \right) - n\hat\mu^2 \\ &= \sum_i \left( x_i^2 - \hat\mu^2 \right) \end{align}$$

which is exactly the "baffling" numerator I encountered in the code.

Improve this answer
$\endgroup$

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.