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 <msu...@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.hard...@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.hard...@wlandres.net> >>> Date: 22-Dec-2013 Time: 09:59:00 >>> This message was sent by XFMail >>> >>> ______________________________________________ >>> R-help@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@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@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.hard...@wlandres.net> Date: 22-Dec-2013 Time: 18:37:15 This message was sent by XFMail ______________________________________________ R-help@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.
