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

[Thomas Lumley]

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