[R] Rounding error in seq(...)

[(Ted Harding)]

On 30-Sep-09 19:32:46, Peter Dalgaard wrote:
> Martin Batholdy wrote:
>> hum,
>> 
>> can you explain that a little more detailed?
>> Perhaps I miss the background knowledge - but it seems just absurd to
>> me.
>> 
>> 0.1+0.1+0.1 is 0.3 - there is no rounding involved, is there?
>> 
>> why is
>> x <- 0.1 + 0.1 +0.1
>> not equal to
>> y <- 0.3
> 
> Remember that this is in BINARY arithmetic. It's really not any
> stranger than the fact that 1/3 + 1/3 != 2/3 in finite accuracy
> decimal arithmetic (0.33333 + 0.33333 = 0.66666 != 0.66667).

Let me perhaps try to spell it out.
Just as in decimal

  1/7 = 0.142857142857142857142857...

with the sequence 142857 recurring for ever (i.e. 1/7 is not
exactly representable in decimal), so

  1/10 = 0.0001100110011001100110011...

in binary, with the sequence 0011 recurring for ever, so 1/10
is not exactly representable in binary (though of course it is
in decimal). When you write (above) "x <- 0.1 + 0.1 +0.1", you
are writing in decimal. R does not do decimal arithmetic (it
could be implemented in a special package perhaps, as it is
in the Unix program 'bc'; but it isn't there), and instead does
its arithmetic internally in binary. Therefore when you enter

  0.1 + 0.1 + 0.1

R will convert each "0.1" into an internal binary representation,
and add these up. On the other hand, when you enter "Y <- 0.3",
R will convert 0.3 into an internal binary representation which
is NOT identical to the result of adding up the three internal
representations of "0.1". These internal representations are to
a fixed number of binary places.

Example in R:
The binary fraction 0.0001100 is (1/2^4 + 1/2^5) (and these are fine
since they are reciprocals of power of 2 and so exactly represented).
The next chunk in the expansion of 1/10 is 0.00011001100 which is
the above, plus 1/16 of itself. And 0.000110011001100 is the above
plus 1/16 plus 1/16^2 of itself, and so on ...

Code:
  M<-1; k<-1/16
  for(i in (1:14)){
    print((1/2^4 + 1/2^5)*M, 17)
    M<-M+k; k<-k/16
  }
# [1] 0.09375
# [1] 0.099609375
# [1] 0.0999755859375
# [1] 0.09999847412109375
# [1] 0.0999999046325684
# [1] 0.0999999940395355
# [1] 0.099999999627471
# [1] 0.0999999999767169
# [1] 0.0999999999985448
# [1] 0.09999999999990905
# [1] 0.0999999999999943
# [1] 0.09999999999999964
# [1] 0.1
# [1] 0.1

(with no further change to 17 decimal places).

Now let's do "0.1 + 0.1 + 0.1" explicitly in binary arithmetic
(to, say, 27 binary places throughout):

  0.000110011001100110011001100 * see below
+ 0.000110011001100110011001100 * see below
-------------------------------
  0.001100110011001100110011000
+ 0.000110011001100110011001100 * see below
-------------------------------
  0.010011001100110011001100100
===============================

which is incorrect, since the last 2 binary digits in the result
should be "10", not "00". The reason is that the next 2 digits in
the truncated representation of 1/10, at the three places marked "*"
above, are "11", and these would get added in if they had been
present. But they are not present, because of truncation, and so
do not get added in. To 27 binary places, the correct representation
of 3/10 (truncated to 27 binary places)is

  0.010011001100110011001100110

Hence (omitting various further details ... )

  3/10 == (1/10+1/10+1/10)
  # [1] FALSE

Because of the cumulative effects of truncation errors like the
above, and rounding errors which can also be introduced along the way,
the method

 (1:9)/10   or 0.1*(1:9)

is better, since each term only acquires error from the specific
arithmetic operation which gives rise to it, and these errors do
not accumulate.

PS: I'm aware of all the grandmothers on the list who do not need
instruction in ovisuction, but this sort of question occurs so
frequently -- and is not spelled out in detail in the docs -- that I
thought it would be useful to actually exhibit an example of the
working for the benefit of people who do wonder just what is going
on then things do not seem to match up as expected!

Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 30-Sep-09                                       Time: 21:53:25
------------------------------ XFMail ------------------------------