R help - Mar 2011 - normal distribution and floating point traps (?): unexpected behavior

dear all,
here's a couple of questions that puzzled me in these last hours:

##### issue 1 qnorm(1-10e-100)!=qnorm(10e-100)

qnorm(1-1e-10) == -qnorm(1e-10)

# turns on to be FALSE. Ok I'm not a computer scientist but,
# but I had a look at the R inferno so I write:

all.equal(qnorm(1-1e-10) , -qnorm(1e-10))

# which turns TRUE, as one would expect, but

all.equal(qnorm(1-1e-100) , -qnorm(1e-100))

# turns FALSE: Indeed

# qnorm(1e-100) is -21.27345, and
# qnorm(1-1e-100) is Inf

# as a consequence there was an infinitive (I would expect a 21.27)
running in my program which was very annoying and hard to detect.

##### issue 2 dnorm(...,log=logical_flag)
##### simple normal likelihood

set.seed(1)
x <- rnorm(100) # sample 100 normal data
mu <- seq(-4,4,length=51) # get a grid for mu
sigma <- seq(.01,4,length=51) # get a grid for sigma
lik <- function(theta,x) sum( dnorm(x,theta[1],theta[2],log=T)) # log
likelihood
Lik <- function(theta,x) prod(dnorm(x,theta[1],theta[2],log=F)) # Likelihood
ms  <- as.matrix(expand.grid(mu,sigma))

l   <- matrix(apply(ms,1,lik,x=x),51)
L   <- matrix(apply(ms,1,Lik,x=x),51)

contour(mu,sigma,L) # plots exactly what I expected
contour(mu,sigma,log(L)) # plots exactly what I expected
contour(mu,sigma,l) # I didn't expect that!!!'
contour(mu,sigma,log(exp(l))) # plots exactly what I expected
> version
               _
platform       x86_64-pc-linux-gnu
arch           x86_64
os             linux-gnu
system         x86_64, linux-gnu
status
major          2
minor          12.2
year           2011
month          02
day            25
svn rev        54585
language       R
version.string R version 2.12.2 (2011-02-25)

under ubuntu distribution

My questions are: is that a bug? Does this behavior depend on my pc
architecture?
What kind of precautions does R's guru suggests to be taken?

Thank you in advance,

PF

ps as usual please I apologize for my fragemnted English

On Wed, Mar 30, 2011 at 10:33 AM, Patrizio Frederic
<frederic.patrizio at gmail.com> wrote:
> dear all,
> here's a couple of questions that puzzled me in these last hours:
>
> ##### issue 1 qnorm(1-10e-100)!=qnorm(10e-100)
>
> qnorm(1-1e-10) == -qnorm(1e-10)
>
> # turns on to be FALSE. Ok I'm not a computer scientist but,
> # but I had a look at the R inferno so I write:
>
> all.equal(qnorm(1-1e-10) , -qnorm(1e-10))
>
> # which turns TRUE, as one would expect, but
>
> all.equal(qnorm(1-1e-100) , -qnorm(1e-100))
>
> # turns FALSE: Indeed
>
> # qnorm(1e-100) is -21.27345, and
> # qnorm(1-1e-100) is Inf
<snip>
> My questions are: is that a bug?
No.  R cannot represent 1-1e-100.  The closest number to 1 that R can
represent is 1-.Machine$double.eps, where .Machine$double.eps is about
1e-16.

FAQ 7.31 is relevant here.
>Does this behavior depend on my pc
> architecture?
In theory yes, though I think all the computers R runs on have the
same value for .Machine$double.eps and certainly none of them can
represent 1-1e-100.
> What kind of precautions does R's guru suggests to be taken?
Use qnorm(1e-100, lower.tail=FALSE) rather than qnorm(1-1e-100)

All the p,q functions have lower.tail= so you can use the more precise
tail of the distribution, and log= so you can get the logarithm of the
result, which can also be useful for increased precision.

    -thomas

-- 
Thomas Lumley
Professor of Biostatistics
University of Auckland

On 29/03/2011 5:33 PM, Patrizio Frederic wrote:
> dear all,
> here's a couple of questions that puzzled me in these last hours:
>
> ##### issue 1 qnorm(1-10e-100)!=qnorm(10e-100)
>
> qnorm(1-1e-10) == -qnorm(1e-10)
>
> # turns on to be FALSE. Ok I'm not a computer scientist but,
> # but I had a look at the R inferno so I write:
>
> all.equal(qnorm(1-1e-10) , -qnorm(1e-10))
>
> # which turns TRUE, as one would expect, but
>
> all.equal(qnorm(1-1e-100) , -qnorm(1e-100))
>
> # turns FALSE: Indeed
>
> # qnorm(1e-100) is -21.27345, and
> # qnorm(1-1e-100) is Inf
Since

1 - 1e-100 == 1

is TRUE, but

1e-100 == 0

is FALSE, this is not a surprise.

Duncan Murdoch