Full_Name: Morten Welinder
Version: 1.5.1
OS: Solaris/Linux
Submission from: (NULL) (192.5.35.38)
The line
return log(1 - p) * (x + 1);
looks like it has problems for p near 1. I would suggest rewriting it to
return log1p (-p) * (x + 1);
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
terra@diku.dk writes:> Full_Name: Morten Welinder > Version: 1.5.1 > OS: Solaris/Linux > Submission from: (NULL) (192.5.35.38) > > > The line > > return log(1 - p) * (x + 1); > > looks like it has problems for p near 1. I would suggest rewriting it to > > return log1p (-p) * (x + 1);I don't think that helps for p near 1. It might do so for p near 0, though. (The whole thing is a bit esoteric since it only becomes relevant when you want the upper tail and on a log scale.) -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Thu, 25 Jul 2002 terra@diku.dk wrote:> Full_Name: Morten Welinder > Version: 1.5.1 > OS: Solaris/Linux > Submission from: (NULL) (192.5.35.38) > > > The line > > return log(1 - p) * (x + 1); > > looks like it has problems for p near 1. I would suggest rewriting it toCould you give us an example of the problems please? Otherwise how will we know that we have solved them?> return log1p (-p) * (x + 1);Please report the bug (see BUGS in the FAQ) and not the solution. It seems to me to be working fine for p as small as 10^(-8) which is already way beyond a feasible application of a geometric (as an exponential would be used instead). -- Brian D. Ripley, ripley@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 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
> return log(1 - p) * (x + 1);> > looks like it has problems for p near 1. I would suggest rewriting it to Could you give us an example of the problems please? Otherwise how will we know that we have solved them? Visual inspection. The problem is that the original code performs an operation, "1-p", with the potential of killing all precision. Often that cannot be helped, but here the function log1p exists as a counterpart to log in order to solve this very problem. > return log1p (-p) * (x + 1); Please report the bug (see BUGS in the FAQ) and not the solution. It seems to me to be working fine for p as small as 10^(-8) which is already way beyond a feasible application of a geometric (as an exponential would be used instead). It is entirely likely that such values are unreasonable, but that is still not a great reason for returning terrible results for something that is at least theoretically sane. The actual point where problems occur is going to depend on your floating point hardware. For me, p=1e-17 makes the old version return zero (fundamentally wrong) and the new version non-zero. Note, that variations of "log (1-p)" exist in other functions, such as pbeta_raw. Morten -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On 26 Jul 2002, Morten Welinder wrote:> > > return log(1 - p) * (x + 1); > > > > looks like it has problems for p near 1. I would suggest rewriting it to > > Could you give us an example of the problems please? Otherwise how > will we know that we have solved them? > > Visual inspection. The problem is that the original code performs an > operation, "1-p", with the potential of killing all precision. OftenBut for p near 1 (you give near 0 below) the distribution is concentrated almost entirely on a single point.> that cannot be helped, but here the function log1p exists as a counterpart > to log in order to solve this very problem. > > > return log1p (-p) * (x + 1); > > Please report the bug (see BUGS in the FAQ) and not the solution. It > seems to me to be working fine for p as small as 10^(-8) which is already > way beyond a feasible application of a geometric (as an exponential would > be used instead). > > It is entirely likely that such values are unreasonable, but that is still > not a great reason for returning terrible results for something that is at > least theoretically sane. The actual point where problems occur is going > to depend on your floating point hardware. For me, p=1e-17 makes the old > version return zero (fundamentally wrong) and the new version non-zero.That's nowhere near one!> Note, that variations of "log (1-p)" exist in other functions, such as > pbeta_raw.I think this amounts to `if I were to do really silly things numerical inaccuracies might catch up with me'. That's true throughout R, but we do have better things to worry about. -- Brian D. Ripley, ripley@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 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._