Thank you for your responses. I presented a minimal example of the issue, but I should have explained that this came up in the context of maximizing a log likelihood function (with optim). I certainly agree that there would be no good reason for a human to evaluate the function pnorm(-x, log.p=TRUE) for large x. However, I would suggest that one may want the function to return a value of -Inf for large x, rather than issue an error or a warning, as I'm guessing your alteration of the C code would do. Returning a value of -Inf would allow an optimization routine (or a pre-optimization-routine grid search) to know that it's going into (or is currently in) the wrong region. I can definitely see the counterargument that someone using an optimization routine should think about the function they are optimizing more carefully, especially given that there could be numerical issues. I think I'm not the right person to judge, both because I don't know the programming issues nor much about numerical optimization, but I wanted to offer a possible argument for returning -Inf in case it hasn't already been considered, and others can evaluate the issue. (In general, my earlier purpose was just to point out what I thought might be an inconsistency in case people wanted to change it. I'm not concerned about it for my own usage at all, since I can just adjust my own programs.)
Thank you, and to all, thank you for all of your work on R. Eric On Fri, 14 May 2010 11:50:09 +0100 (BST), Prof Brian Ripley <rip...@stats.ox.ac.uk> wrote: > 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.hard...@manchester.ac.uk wrote: > >> On 13-May-10 20:04:50, efree...@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.hard...@manchester.ac.uk> >> Fax-to-email: +44 (0)870 094 0861 >> Date: 14-May-10 Time: 00:01:27 >> ------------------------------ XFMail ------------------------------ >> >> ______________________________________________ >> R-devel@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-devel >> ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
