[R] normal distribution and floating point traps (?): unexpected behavior

Patrizio Frederic frederic.patrizio at gmail.com
Tue Mar 29 23:33:42 CEST 2011


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



More information about the R-help mailing list