>>>>> "daheiser" == daheiser <daheiser@gvn.net>
>>>>> on Tue, 6 Apr 2004 04:24:35 +0200 (CEST) writes:
daheiser> It is an error in the algorithm.
"it" being the behavior reported in bug report PR#1228 ---
[too bad you didn't use the whole string "PR#1228" in your
subject;
if you had, no new report would have been created, and things
would have gone to the proper place in the R-bugs repository ...]
namely the fact that
var(rep(1e30, n))
does not always (for all n) give 0. With the following function
tst.var <- function(x,nset= 2:100) {
for(n in nset) {
v <- var(rep(x,n))
if(v != 0) cat(sprintf("%4d: %20.12g\n", n,v))
}
}
Actually, on my current desktop (AMD Athlon,Linux, R 1.9.0beta)
the computations seem to be more often exact than they used to:
Even tst.var(1e30, 2:5000) doesn't produce any output nor do a few other
'x'
arguments I tried --- but the fundamental criticism is still
correct, and even on the Athlon, it's easy to find cases where
tst.var() *does* produce output.
daheiser> It is an excellent example of how errors and
daheiser> faults occur when the programmer follows the
daheiser> mathematical formula exactly. Welford's algorithm
daheiser> does not produce this error. It gives correct
daheiser> standard deviation and variance values.
Actually, after reading
@article{ChaTL97,
author = {Tony F. Chan and John Gregg Lewis},
title = {Computing standard deviations: accuracy},
journal = {Commun. ACM},
volume = 22,
number = 9,
year = 1979,
issn = {0001-0782},
pages = {526--531},
doi = {http://doi.acm.org/10.1145/359146.359152},
publisher = {ACM Press},
}
@article{WesD97,
author = {D. H. D. West},
title = {Updating mean and variance estimates: an improved method},
journal = {Commun. ACM},
volume = 22,
number = 9,
year = 1979,
issn = {0001-0782},
pages = {532--535},
doi = {http://doi.acm.org/10.1145/359146.359153},
publisher = {ACM Press},
}
I'd conclude (from page 531) that Welford's algorithm is a bit
less accurate than the (very similar) "West" version,
and we (the R developers) should rather implement the latter.
Maybe for R 1.9.1 -- or even later {there are even more
important numerical accuracy problems I know about!}
Martin Maechler <maechler@stat.math.ethz.ch>
http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum LEO C16 Leonhardstr. 27
ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND
phone: x-41-1-632-3408 fax: ...-1228 <><
