[Rd] Computing means, variances and sums

[Prof Brian Ripley]

On Sun, 19 Feb 2006, hadley wickham wrote:

>> p.s.  If my computations are correct, 0.2 = 0*/2 + 0/4 + 1/8 + 1/16 +
>> 0/32 + 0/64 + 1/128 + 1/256 + 0/512 + 0/1024 + 1/2048 + 1/4096 + ... =
>> 0.3333333333333h.  Perhaps someone can extend this to an FAQ to help
>> explain finite precision arithmetic and rounding issues.
>
> This is drifting a bit off topic, but the other day I discovered this
> rather nice illustration of the perils of finite precision arithmetic
> while creating a contrast matrix:
>
>> n <- 13
>> a <- matrix(-1/n, ncol=n, nrow=n) + diag(n)
>> rowSums(a)
> [1]  2.775558e-16  2.775558e-16  5.551115e-17  5.551115e-17  5.551115e-17
> [6]  5.551115e-17  0.000000e+00 -5.551115e-17  0.000000e+00  5.551115e-17
> [11]  1.110223e-16  1.665335e-16  2.220446e-16
>
> Not only do most of the rows not sum to 0, they do not even sum to the
> same number!  It is hard to remember the familiar rules of arithmetic
> do not always apply.

I think you will find this example does give all 0's in R-devel, even 
on platforms like Sparc.  But users do need to remember that computer 
arithmetic is inexact except in rather narrowly delimited cases.

-- 
Brian D. Ripley,                  ripley at 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 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595