I've been following this thread with interest. A nice
collection of things to watch out for, if you don't
want the small arithmetic errors due to finite-length
digital representations of fractions to cause trouble!
However, as well as these small discrepancies, major
malfunctions can also result.
Back on Dec 22, 2013, I posted a Christmas Greetings
message to R-help:
Season's Greetings (and great news ... )!
which starts:
Greetings All!
With the Festive Season fast approaching, I bring you joy
with the news (which you will surely wish to celebrate)
that R cannot do arithmetic!
Usually, this is manifest in a trivial way when users report
puzzlement that, for instance,
sqrt(pi)^2 == pi
# [1] FALSE
which is the result of a (generally trivial) rounding or
truncation error:
sqrt(pi)^2 - pi
# [1] -4.440892e-16
But for some very simple calculations R goes off its head.
And the example given is:
Consider a sequence generated by the recurrence relation
x[n+1] = 2*x[n] if 0 <= x[n] <= 1/2
x[n+1] = 2*(1 - x[n]) if 1/2 < x[n] <= 1
(for 0 <= x[n] <= 1).
This has equilibrium points (x[n+1] = x[n]) at x[n] = 0
and at x[n] = 2/3:
2/3 -> 2*(1 - 2/3) = 2/3
It also has periodic points, e.g.
2/5 -> 4/5 -> 2/5 (period 2)
2/9 -> 4/9 -> 8/9 -> 2/9 (period 3)
The recurrence relation can be implemented as the R function
nextx <- function(x){
if( (0<=x)&(x<=1/2) ) {x <- 2*x} else {x <- 2*(1 - x)}
}
Now have a look at what happens when we start at the equilibrium
point x = 2/3:
N <- 1 ; x <- 2/3
while(x > 0){
cat(sprintf("%i: %.9f\n",N,x))
x <- nextx(x) ; N <- N+1
}
cat(sprintf("%i: %.9f\n",N,x))
For a while [run it and see!], this looks as though it's doing what
the arithmetic would lead us to expect: the first 24 results will all
be printed as 0.666666667, which looks fine as 2/3 to 9 places.
But then the "little errors" start to creep in:
N=25: 0.666666666
N=28: 0.666666672
N=46: 0.667968750
N=47: 0.664062500
N=48: 0.671875000
N=49: 0.656250000
N=50: 0.687500000
N=51: 0.625000000
N=52: 0.750000000
N=53: 0.500000000
N=54: 1.000000000
N=55: 0.000000000
What is happening is that, each time R multiplies by 2, the binary
representation is shifted up by one and a zero bit is introduced
at the bottom end.
At N=53, the first binary bit of 'x' is 1, and all the rest are 0,
so now 'x' is exactly 0.5 = 1/2, hence the final two are also exact
results; 53 is the Machine$double.digits = 53 binary places.
So this normally "almost" trivial feature can, for such a simple
calculation, lead to chaos or catastrophe (in the literal technical
sense).
For more detail, including an extension of the above, look at the
original posting in the R-help archives for Dec 22, 2013:
From: (Ted Harding) <ted.hard...@wlandres.net>
Subject: [R] Season's Greetings (and great news ... )!
Date: Sun, 22 Dec 2013 09:59:16 -0000 (GMT)
(Apologies, but I couldn't track down the URL for this posting
in the R-help archives; there were a few follow-ups).
I gave this as an example to show that the results of the "little"
arithmetic errors (such as have recently been discussed from many
aspects) can, in certain contexts, destroy a computation.
So be careful to consider what can happen in the particular
context you are working with.
There are ways to dodge the issue -- such as using the R interface
to the 'bc' calculator, which computes arithmetic expressions in
a way which is quite different from the fixed-finite-length binary
representation and algorithms used, not only by R, but also by many
other numerical computation software suites
Best wishes to all,
Ted.
-------------------------------------------------
E-Mail: (Ted Harding) <ted.hard...@wlandres.net>
Date: 21-Apr-2017 Time: 22:03:15
This message was sent by XFMail
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
