[Rd] bug with var(rep(1e30, 3)) (PR#1228)
ripley@stats.ox.ac.uk
ripley@stats.ox.ac.uk
Wed, 26 Dec 2001 16:52:57 +0100 (MET)
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.
On 26 Dec 2001, Peter Dalgaard BSA wrote:
> paradis@isem.univ-montp2.fr writes:
>
> > There seems to be a bug with var() when the argument is a vector with
> > exactly three values of 1e30 (or close to this value). This does not happen
> > with twice, four (or more) times this value, or another value.
> >
> >
> > > var(rep(1e30, 3))
> > [1] 2.971056e+28
> > > var(rep(1.2e30, 3))
> > [1] 2.971056e+28
> > > var(rep(0.9e30, 3))
> > [1] 2.971056e+28
> > > var(rep(0.8e30, 3))
> > [1] 0
> > > var(rep(1e29, 3))
> > [1] 0
> > > var(rep(1e30, 2))
> > [1] 0
> > > var(rep(1e30, 4))
> > [1] 0
> > > var(rep(1e31, 3))
> > [1] 0
> >
> >
> > The bug is repeatable, and I got the same results with R 1.3.1 (binary from
> > CRAN) and R 1.4.0 (binary from BD Ripley's site) both on Windows NT, and R
> > 1.4.0 on Solaris (compiled from sources).
>
> Not a bug, roundoff:
>
> > (1e30+1e30+1e30)/3-1e30
> [1] 1.407375e+14
> > 3*((1e30+1e30+1e30)/3-1e30)^2/2
> [1] 2.971056e+28
>
> Relative errors of the order 1e-16 are completely expectable in
> floating point arithmetic.
>
> For illustration, look at it in octal, after removing a bunch of
> trailing zeroes:
>
> > 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
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>
--
Brian D. Ripley, ripley@stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272860 (secr)
Oxford OX1 3TG, UK Fax: +44 1865 272595
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._