[Rd] bug with var(rep(1e30, 3)) (PR#1228)

[p.dalgaard@biostat.ku.dk]

Prof Brian Ripley <ripley@stats.ox.ac.uk> writes:

> But there are better algorithms that are not subject to this sort
> of rounding error.
> 
> I reckon that computing a variance this way is not the quality of
> algorithm one should expect from R.

Hmm. Which algorithms are there? Surely almost any kind of provisional
means subtraction could avoid the the extreme case where the mean of a
bunch of identical numbers is different from all of them, but will it
also improve precision in the general case?

Essentially our current algorithm in cov.c is

xm<- sum(x)/n
v <- 1/(n-1)*sum((x-xm)^2)

> > > x<-1e30/(8^16)
> > > structure(x, class="octmode")
> > [1] "144762623464043165"
> > > structure(x+x+x, class="octmode")
> > [1] "456730272634151540"
> >
> > If you look carefully, you will see that a bit is lost in the process
> > since 5+5+5 is 17 octal.


-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._