[Rd] results of pnorm as either NaN or Inf

[Prof Brian Ripley]

The answer to pnorm(-x, log.p=TRUE) is of course about -0.5*x*x for 
large x.  So trying to compute it for -1e108 is a bit silly, both on 
your part and the C code used.  I've altered the C code not to attempt 
it for x > 1e170, which excludes the area where underflow occurs.

On Fri, 14 May 2010, Ted.Harding at manchester.ac.uk wrote:

> On 13-May-10 20:04:50, efreeman at berkeley.edu wrote:
>> I stumbled across this and I am wondering if this is unexpected
>> behavior or if I am missing something.
>>
>>> pnorm(-1.0e+307, log.p=TRUE)
>> [1] -Inf
>>> pnorm(-1.0e+308, log.p=TRUE)
>> [1] NaN
>> Warning message:
>> In pnorm(q, mean, sd, lower.tail, log.p) : NaNs produced
>>> pnorm(-1.0e+309, log.p=TRUE)
>> [1] -Inf
>>
>> I don't know C and am not that skilled with R, so it would be hard
>> for me to look into the code for pnorm. If I'm not just missing
>> something, I thought it may be of interest.
>>
>> Details:
>> I am using Mac OS X 10.5.8. I installed a precompiled binary version.
>> Here is the output from sessionInfo(), requested in the posting guide:
>> R version 2.11.0 (2010-04-22)
>> i386-apple-darwin9.8.0
>>
>> locale:
>> [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
>>
>> attached base packages:
>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>
>> loaded via a namespace (and not attached):
>> [1] tools_2.11.0
>>
>> Thank you very much,
>>
>> Eric Freeman
>> UC Berkeley
>
> This is probably platform-independent. I get the same results with
> R on Linux. More to the point:
>
> You are clearly "pushing the envelope" here. First, have a look
> at what R makes of your inputs to pnorm():
>
>  -1.0e+307
>  # [1] -1e+307
>  -1.0e+308
>  # [1] -1e+308
>  -1.0e+309
>  # [1] -Inf
>
>
> So, somewhere beteen -1e+308 and -1.0e+309. the envelope burst!
> Given -1.0e+309, R returns -Inf (i.e. R can no longer represent
> this internally as a finite number).
>
> Now look at
>
>  pnorm(-Inf,log.p=TRUE)
>  # [1] -Inf
>
> So, R knows how to give the correct answer (an exact 0, or -Inf
> on the log scale) if you feed pnorm() with -Inf. So you're OK
> with -1e+N where N >= 309.
>
> For smaller powers, e.g. -1e+(200:306), these give pnorm() much
> less than -1.0e+309, and presumably R's algorithm (which I haven't
> studied either ... ) returns 0 for pnorm(), as it should to the
> available internal accuracy.
>
> So, up to pnorm(-1.0e+307, log.p=TRUE) = -Inf. All is as it should be.
> However, at -1e+308, "the envelope is about to burst", and something
> may occur within the algorithm which results in a NaN.
>
> So there is nothing anomalous about your results except at -1e+308,
> which is where R is at a critical point.
>
> That's how I see it, anway!
> Ted.
>
> --------------------------------------------------------------------
> E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
> Fax-to-email: +44 (0)870 094 0861
> Date: 14-May-10                                       Time: 00:01:27
> ------------------------------ XFMail ------------------------------
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595