Re: [R] Season's Greetings (and great news ... )!
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.7
2: 0.7
3: 0.7
4: 0.7
5: 0.7
...
but, later on,
24: 0.7
25: 0.6
26: 0.8
27: 0.4
28: 0.66672
...
46: 0.667968750
47: 0.664062500
48: 0.671875000
49: 0.65625
50: 0.68750
51: 0.62500
52: 0.75000
53: 0.5
54: 1.0
55: 0.0
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.000 - 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.101 x[5] - 2*(1 - 0.101) = 2*0.011 -
0.110 x[6] - 2*(1 - 0.110) = 2*0.010 -
0.100 x[7] - 2*0.100 = 1.000 -
1.000 x[8] - 2*(1 - 1.000) = 2*0 -
0.000 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(x0){
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
Re: [R] Season's Greetings (and great news ... )!
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.7
2: 0.7
3: 0.7
4: 0.7
5: 0.7
...
but, later on,
24: 0.7
25: 0.6
26: 0.8
27: 0.4
28: 0.66672
...
46: 0.667968750
47: 0.664062500
48: 0.671875000
49: 0.65625
50: 0.68750
51: 0.62500
52: 0.75000
53: 0.5
54: 1.0
55: 0.0
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.000 - 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.101 x[5] - 2*(1 - 0.101) = 2*0.011 -
0.110 x[6] - 2*(1 - 0.110) = 2*0.010 -
0.100 x[7] - 2*0.100 = 1.000 -
1.000 x[8] - 2*(1 - 1.000) = 2*0 -
0.000 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(x0){
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)
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 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.7
2: 0.7
3: 0.7
4: 0.7
5: 0.7
...
but, later on,
24: 0.7
25: 0.6
26: 0.8
27: 0.4
28: 0.66672
...
46: 0.667968750
47: 0.664062500
48: 0.671875000
49: 0.65625
50: 0.68750
51: 0.62500
52: 0.75000
53: 0.5
54: 1.0
55: 0.0
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.000 - 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.101 x[5] - 2*(1 - 0.101) = 2*0.011 -
0.110 x[6] - 2*(1 - 0.110) = 2*0.010 -
0.100 x[7] - 2*0.100 = 1.000 -
1.000 x[8] - 2*(1 - 1.000) = 2*0 -
0.000 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 ;
Re: [R] Season's Greetings (and great news ... )!
(or whatever the keyboard analogue of that may be)
Hands clasped? Fingers interlaced?
John Kane
Kingston ON Canada
-Original Message-
From: ted.hard...@wlandres.net
Sent: Sun, 22 Dec 2013 18:37:18 - (GMT)
To: r-help@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 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.7
2: 0.7
3: 0.7
4: 0.7
5: 0.7
...
but, later on,
24: 0.7
25: 0.6
26: 0.8
27: 0.4
28: 0.66672
...
46: 0.667968750
47: 0.664062500
48: 0.671875000
49: 0.65625
50: 0.68750
51: 0.62500
52: 0.75000
53: 0.5
54: 1.0
55: 0.0
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.000 - 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.101 x[5] - 2*(1 - 0.101) = 2*0.011 -
0.110 x[6] - 2*(1 - 0.110) = 2*0.010 -
0.100 x[7] - 2*0.100 = 1.000 -
1.000 x[8] - 2*(1 - 1.000) = 2*0 -
0.000 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
