[R] Season's Greetings (and great news ... )!

John Kane jrkrideau at inbox.com
Sun Dec 22 19:41:11 CET 2013


"(or whatever the keyboard analogue of that may be) "

Hands clasped?  Fingers interlaced?

John Kane
Kingston ON Canada


> -----Original Message-----
> From: ted.harding at wlandres.net
> Sent: Sun, 22 Dec 2013 18:37:18 -0000 (GMT)
> To: r-help at r-project.org
> Subject: Re: [R] Season's Greetings (and great news ... )!
> 
> Thanks for the comments, Bert and Mehmet! It is of course a serious
> and interesting area to explore (and I'm aware of the chaos context;
> I initially got into this areas year ago when I was exploring the
> possibilities for chaos in fish population dynamics -- and they're
> certainly there)!
> 
> But, before anyone takes my posting *too* seriously, let me say that
> it was written tongue-in-cheek (or whatever the keyboard analogue of
> that may be). I'm certainly not "blaming R".
> 
> Have fun anyway!
> Ted.
> 
> On 22-Dec-2013 17:35:56 Bert Gunter wrote:
>> Yes.
>> 
>> See also Feigenbaum's constant and chaos theory for the general context.
>> 
>> Cheers,
>> Bert
>> 
>> On Sun, Dec 22, 2013 at 8:54 AM, Suzen, Mehmet <msuzen at gmail.com> wrote:
>>> I wouldn't blame R for floating-point arithmetic and our personal
>>> feeling of what 'zero' should be.
>>> 
>>>> options(digits=20)
>>>> pi
>>> [1] 3.141592653589793116
>>>> sqrt(pi)^2
>>> [1] 3.1415926535897926719
>>>> (pi - sqrt(pi)^2) < 1e-15
>>> [1] TRUE
>>> 
>>> There was a similar post before, for example see:
>>> http://r.789695.n4.nabble.com/Why-does-sin-pi-not-return-0-td4676963.html
>>> 
>>> There is an example by Martin Maechler (author of Rmpfr) on how to use
>>> arbitrary precision
>>> with your arithmetic.
>>> 
>>> On 22 December 2013 10:59, Ted Harding <Ted.Harding at wlandres.net>
>>> wrote:
>>>> 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.
>>>> 
>>>> I had originally posted this example some years ago, but I
>>>> have since generalised it, and the generalisation is even
>>>> more entertaining than the original.
>>>> 
>>>> The Original:
>>>> 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))
>>>> 
>>>> Run that, and you will see that successive values of x collapse
>>>> towards zero. Things look fine to start with:
>>>> 
>>>>   1: 0.666666667
>>>>   2: 0.666666667
>>>>   3: 0.666666667
>>>>   4: 0.666666667
>>>>   5: 0.666666667
>>>>   ...
>>>> 
>>>> but, later on,
>>>> 
>>>>   24: 0.666666667
>>>>   25: 0.666666666
>>>>   26: 0.666666668
>>>>   27: 0.666666664
>>>>   28: 0.666666672
>>>>   ...
>>>> 
>>>>   46: 0.667968750
>>>>   47: 0.664062500
>>>>   48: 0.671875000
>>>>   49: 0.656250000
>>>>   50: 0.687500000
>>>>   51: 0.625000000
>>>>   52: 0.750000000
>>>>   53: 0.500000000
>>>>   54: 1.000000000
>>>>   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. To illustrate this, do the calculation in
>>>> 7-bit arithmetic where 2/3 = 0.1010101, so:
>>>> 
>>>> 0.1010101  x[1], >1/2 so subtract from 1 = 1.0000000 -> 0.0101011,
>>>> and then multiply by 2 to get x[2] = 0.1010110. Hence
>>>> 
>>>> 0.1010101  x[1] -> 2*(1 - 0.1010101) = 2*0.0101011 ->
>>>> 0.1010110  x[2] -> 2*(1 - 0.1010110) = 2*0.0101010 ->
>>>> 0.1010100  x[3] -> 2*(1 - 0.1010100) = 2*0.0101100 ->
>>>> 0.1011000  x[4] -> 2*(1 - 0.1011000) = 2*0.0101000 ->
>>>> 0.1010000  x[5] -> 2*(1 - 0.1010000) = 2*0.0110000 ->
>>>> 0.1100000  x[6] -> 2*(1 - 0.1100000) = 2*0.0100000 ->
>>>> 0.1000000  x[7] -> 2*0.1000000 = 1.0000000 ->
>>>> 1.0000000  x[8] -> 2*(1 - 1.0000000) = 2*0 ->
>>>> 0.0000000  x[9] and the end of the line.
>>>> 
>>>> The final index of x[i] is i=9, 2 more than the number of binary
>>>> places (7) in this arithmetic, since 8 successive zeros have to
>>>> be introduced. It is the same with the real R calculation since
>>>> this is working to .Machine$double.digits = 53 binary places;
>>>> it just takes longer (we reach 0 at x[55])! The above collapse
>>>> to 0 occurs for any starting value in this simple example (except
>>>> for multiples of 1/(2^k), when it works properly).
>>>> 
>>>> Generalisation:
>>>> This is basically the same, except that everything is multiplied
>>>> by a scale factor S, so instead of being on the interval [0,1].
>>>> it is on [0,S], and
>>>> 
>>>>   x[n+1] = 2*x[n] if 0 <= x[n] <= S/2
>>>>   x[n+1] = 2*(S - x[n]) if S/2 < x[n] <= S
>>>> (for 0 <= x[n] <= S).
>>>> 
>>>> Again, x[n] = 2*S/3 is an equilibrium point. 2*S/3 > S/2, so
>>>> 
>>>>   x[n] -> 2*(S - 2*S/3) = 2*(S/3) = 2*S/3
>>>> 
>>>> Functions to implement this:
>>>> 
>>>>   nxtS <- function(x,S){
>>>>     if((x >= 0)&(x <= S/2)){ x<- 2*x } else {x <- 2*(S-x)}
>>>>   }
>>>> 
>>>>   S <- 6 ##  Or some other value of S
>>>>   Nits <- 100
>>>>   x <- 2*S/3
>>>>   N <- 1 ; print(c(N,x))
>>>>   while(x>0){
>>>>   if(N > Nits) break   ### to stop infinite looping
>>>>   N <- (N+1) ; x <- nxtS(x,S)
>>>>   print(c(N,x))
>>>> }
>>>> 
>>>> The behaviour of the sequence now depends on the value of S.
>>>> 
>>>> If S is a multiple of 3, then with x[1] = 2*S/3 the equilibrium
>>>> is immediately attained and x[n] = 2*S/3 forever after, since
>>>> R is now calculating with integers. E.g. try the above with S<-6
>>>> That is what arithmetic ought to be like! But for S not a multiple
>>>> of 3 one can get the impression that R is on some sort of drug!
>>>> 
>>>> For other values of S (but not all) we observe the same collapse
>>>> to x=0 as before, and again it takes 54 steps (ending with x[55]).
>>>> Try e.g. S <- 16
>>>> 
>>>> For some values of S, however, the iteration ends up in a periodic
>>>> loop.
>>>> 
>>>> For example, with S<-7, at x[52] we get x[52]=4, x[53]=6, x[54]=2,
>>>> and then 4 6 2 4 6 2 4 6 2 ... forever (or until Nits cuts in),
>>>> so period = 3.
>>>> 
>>>> For S<-11, x[52]=8 then 6 then 10 then 2 then 4 then 8 6 10 2 4 ...
>>>> so period = 5.
>>>> 
>>>> For S<-13, x[51]=4 then 8 10 6 12 2 4 8 10 6 12 2 4 8 ...
>>>> so period = 6.
>>>> 
>>>> For S<-19, x[51]=12 then 14 10 18 2 4 8 16 6 12 ...
>>>> so period = 9.
>>>> 
>>>> And so on ...
>>>> 
>>>> So, one sniff of something like S<-19, and R is off its head!
>>>> 
>>>> All it has to do is multiply by 2 -- and it gets it cumulatively
>>>> wrong!
>>>> R just doesn't add up ...
>>>> 
>>>> Season's Greetings to all -- and may your calculations always
>>>> be accurate -- to within machine precision ...
>>>> 
>>>> Ted.
>>>> 
>>>> -------------------------------------------------
>>>> E-Mail: (Ted Harding) <Ted.Harding at wlandres.net>
>>>> Date: 22-Dec-2013  Time: 09:59:00
>>>> This message was sent by XFMail
>>>> 
>>>> ______________________________________________
>>>> R-help at r-project.org mailing list
>>>> 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.
>>> 
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> 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.
>> 
>> 
>> 
>> --
>> 
>> Bert Gunter
>> Genentech Nonclinical Biostatistics
>> 
>> (650) 467-7374
>> 
>> ______________________________________________
>> R-help at r-project.org mailing list
>> 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.
> 
> -------------------------------------------------
> E-Mail: (Ted Harding) <Ted.Harding at wlandres.net>
> Date: 22-Dec-2013  Time: 18:37:15
> This message was sent by XFMail
> 
> ______________________________________________
> R-help at r-project.org mailing list
> 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.

____________________________________________________________
FREE 3D MARINE AQUARIUM SCREENSAVER - Watch dolphins, sharks & orcas on your desktop!



More information about the R-help mailing list