src/nmath/pgeom.c (PR#1834)
ripley@stats.ox.ac.uk
ripley@stats.ox.ac.uk
Fri, 26 Jul 2002 17:16:54 +0200 (MET DST)
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. Often
But 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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._