src/nmath/pgeom.c (PR#1834)

[terra@diku.dk]

   >     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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._