I stumbled across a paper by Welford (1962), where he proclaims a method that should compute the standard deviation numerically more robust than the naive algorithms (http://www.jstor.org/stable/1266577).
Here, "numerically robust" means that round-off errors are reduced.
He gives a recurrence for the sum of squares $S_n = \sum_{i=1}^n (x_i – \mu_n)^2 = S_{n-1} + \frac{n-1}{n} (x_n – \mu_{n-1})^2 $ , $\mu$ being, of course, the mean.
That way, the standard deviation can be computed iteratively in a single pass.
As far as I understand, he claims that his iterative recurrence formula is numerically more robust, because all terms in it are of the same order (provided the input data are all of the same order).
This is what I don't understand. It seems to me that, as $n$ gets incremented, $S_n$ becomes larger and larger, so more and more significant digits from the second term are lost, aren't they?
Googling a bit further, I have found a paper by
Youngs & Cramer, 1971, who looked at a number of methods, including Welford's, of computing the sum of products / standard deviation more robustly (http://www.jstor.org/stable/1267176).
Conducting a number of experiments, they found that Welford's method does not provide any benefits. So that seems to confirm my doubts about Welford's method.
Now, Youngs & Cramer propose another method, which computes $S_n = S_{n-1} + \frac{n-1}{n}(nx_n – s_n)^2 $, where $s_n = s_{n-1} + x_n$.
Empirically, they found their method to be superior.
Again, I don't understand why this should be the case: isn't there some catastrophic cancellation going on in $(nx_n – s_n)$ ? Don't the terms $S_n$ and the fraction differ in their magnitude more and more, so that more and more digits of the fraction term get rounded off?
I would by most grateful if somebody could shed some light on these questions.
In addition, I'd like to know which is the best method (in terms of roundoff errors) to compute these sums. Surprisingly, I haven't found anything about this in Numerical Recipes (or I overlooked it).
Thank you very much in advance.
Gabriel.
Best Answer
I found a discussion of this exact problem in Higham, Accuracy and stability of numerical algorithms, Section 1.9.
The author suggests an alternative algorithm and claims that it is numerically stable (in the mixed backward-forward sense); proofs for the precise accuracy bounds of the formulas are left as exercises.