From murdoch at stats.uwo.ca  Sun Feb 19 22:58:23 2006
From: murdoch at stats.uwo.ca (Duncan Murdoch)
Date: Sun, 19 Feb 2006 16:58:23 -0500
Subject: [Rd] Computing means, variances and sums
In-Reply-To: <Pine.LNX.4.64.0602192014300.28769@gannet.stats.ox.ac.uk>
References: <Pine.LNX.4.64.0602190957290.4741@gannet.stats.ox.ac.uk>	<43F8C0D1.6000907@pdf.com>	<f8e6ff050602191142x6a7c1359v768980d0f8842bf2@mail.gmail.com>
	<Pine.LNX.4.64.0602192014300.28769@gannet.stats.ox.ac.uk>
Message-ID: <43F8E9FF.9080304@stats.uwo.ca>

On 2/19/2006 3:18 PM, Prof Brian Ripley wrote:
> On Sun, 19 Feb 2006, hadley wickham wrote:
> 
>>> p.s.  If my computations are correct, 0.2 = 0*/2 + 0/4 + 1/8 + 1/16 +
>>> 0/32 + 0/64 + 1/128 + 1/256 + 0/512 + 0/1024 + 1/2048 + 1/4096 + ... =
>>> 0.3333333333333h.  Perhaps someone can extend this to an FAQ to help
>>> explain finite precision arithmetic and rounding issues.
>> This is drifting a bit off topic, but the other day I discovered this
>> rather nice illustration of the perils of finite precision arithmetic
>> while creating a contrast matrix:
>>
>>> n <- 13
>>> a <- matrix(-1/n, ncol=n, nrow=n) + diag(n)
>>> rowSums(a)
>> [1]  2.775558e-16  2.775558e-16  5.551115e-17  5.551115e-17  5.551115e-17
>> [6]  5.551115e-17  0.000000e+00 -5.551115e-17  0.000000e+00  5.551115e-17
>> [11]  1.110223e-16  1.665335e-16  2.220446e-16
>>
>> Not only do most of the rows not sum to 0, they do not even sum to the
>> same number!  It is hard to remember the familiar rules of arithmetic
>> do not always apply.
> 
> I think you will find this example does give all 0's in R-devel, even 
> on platforms like Sparc.  

Only until the fpu precision gets changed:

 > n <- 13
 > a <- matrix(-1/n, ncol=n, nrow=n) + diag(n)
 > rowSums(a)
  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0
 > RSiteSearch('junk')
A search query has been submitted to http://search.r-project.org
The results page should open in your browser shortly

 > n <- 13
 > a <- matrix(-1/n, ncol=n, nrow=n) + diag(n)
 > rowSums(a)
  [1]  2.775558e-16  2.775558e-16  5.551115e-17  5.551115e-17  5.551115e-17
  [6]  5.551115e-17  0.000000e+00 -5.551115e-17  0.000000e+00  5.551115e-17
[11]  1.110223e-16  1.665335e-16  2.220446e-16

We still need to protect against these changes.  I'll put something 
together, unless you're already working on it.

The approach I'm thinking of is to define a macro to be called in risky 
situations.  On platforms where this isn't an issue, the macro would be 
null; on Windows, it would reset the fpu to full precision.

For example, RSiteSearch causes damage in the ShellExecute call in 
do_shellexec called from browseURL, so I'd add protection there.  I 
think we should also add detection code somewhere in the evaluation loop 
to help in diagnosing these problems.

> But users do need to remember that computer 
> arithmetic is inexact except in rather narrowly delimited cases.

Yes, that too.

Duncan Murdoch



[Rd] Computing means, variances and sums
Duncan Murdoch