On 04/21/2017 02:03 PM, (Ted Harding) wrote:
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).
The only surprise is to end up at the other equilibrium point (0),
not that the iterations didn't stabilize around or converge to
equilibrium point 2/3. That's because:
1) You didn't start *exactly* at equilibrium point 2/3
2) The iterating process is not supposed to converge but to diverge
from the equilibrium point. It's its mathematical nature and has
nothing to do with rounding errors.
So if you don't start exactly at 2/3, the iterations will take you
further apart, even on a hypothetical computer that uses exact
representation of any arbitrary real number (i.e. no rounding errors).
With such an instable algorithm, all bets are off: you can end up
to the other equilibrium point or not end up anywhere (chaos) or
hit a periodic point and go in cycles from there.
But hey, shouldn't people examine the mathematical properties of
their iterative numerical algorithms before implementing them? Like
convergence, stability etc...
H.
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://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mailman_listinfo_r-2Dhelp&d=DwICAg&c=eRAMFD45gAfqt84VtBcfhQ&r=BK7q3XeAvimeWdGbWY_wJYbW0WYiZvSXAJJKaaPhzWA&m=ILkgA-W4NkL_PpalU1VtS9iO60moznxYVBGfU5QbyNw&s=puQK7QZvG7lZuOYvo9D6UANEgWylkApF-xIvpLMGVhQ&e=
PLEASE do read the posting guide
https://urldefense.proofpoint.com/v2/url?u=http-3A__www.R-2Dproject.org_posting-2Dguide.html&d=DwICAg&c=eRAMFD45gAfqt84VtBcfhQ&r=BK7q3XeAvimeWdGbWY_wJYbW0WYiZvSXAJJKaaPhzWA&m=ILkgA-W4NkL_PpalU1VtS9iO60moznxYVBGfU5QbyNw&s=zohn8q-ufE7UhlQDQGc3KP-lsbM4O62SzwBXQ4S7I0g&e=
and provide commented, minimal, self-contained, reproducible code.
--
Hervé Pagès
Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M1-B514
P.O. Box 19024
Seattle, WA 98109-1024
E-mail: hpa...@fredhutch.org
Phone: (206) 667-5791
Fax: (206) 667-1319
______________________________________________
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.
