R devel - Aug 2000 - Re: [R] too large alpha or beta in dbeta ? (PR#643)

>>>>> "MM" == Martin Maechler
<maechler@stat.math.ethz.ch> writes:
>>>>> "TL" == Thomas Lumley
<thomas@biostat.washington.edu> writes:
    TL> On Thu, 24 Aug 2000, Troels Ring wrote:
    >>> Dear friends.
    >>> 
    >>> Is this as expected ? Is alpha and beta too large simply ?
    >>> 
    >>> > dbeta(.1,534,646)
    >>> [1] NaN
    >>> Warning message:
    >>> NaNs produced in: dbeta(x, shape1, shape2, log)

    TL> well, it should work, but the correct answer is effectively zero.
    TL> pbeta(.1,534,646) gives 3.6e-213

    TL> Perhaps more worrying is 
    >>> dbeta(.25,534,646)
    TL> [1] Inf
    MM> yes.
    MM> and I see that it is one case where the log-density seems to be
alright:

    >> dbeta(.25,534,646,log=TRUE)
    MM> [1] -109.939612

    MM> and also for the NaN case :

    >> dbeta(.1,534,646, log=TRUE)
    MM> [1] -480.725168

    MM> A workaround is  using   exp( log-density  ), i.e.

    MM> exp(dbeta(x,a,b, log = TRUE)) :

    MM> Look at

    >> plot(function(x)    dbeta(x, 534,646, log = TRUE), n = 1001)
    MM> or
    >> plot(function(x)exp(dbeta(x, 534,646, log = TRUE)), n = 1001)

    MM> --
    MM> I'll have a look.

fixed for R-devel (1.2 unstable),

and the patch is

--- src/nmath/dbeta.c	2000/03/03 17:18:30	1.7
+++ src/nmath/dbeta.c	2000/08/25 16:16:35	1.8
@@ -30,7 +30,15 @@
 #ifdef IEEE_754
     /* NaNs propagated correctly */
     if (ISNAN(x) || ISNAN(a) || ISNAN(b)) return x + a + b;
+
+# define xmax 171.61447887182298/* (fixme) -->> ./gammalims.c */
+
+#else
+    static double xmax = 0; double xmin;
+    if (xmax == 0)
+	gammalims(&xmin, &xmax);
 #endif
+
     if (a <= 0 || b <= 0) ML_ERR_return_NAN;
 
     if (x < 0 || x > 1)
@@ -35,10 +43,13 @@
 
     if (x < 0 || x > 1)
 	return R_D__0;
+
+#define R_LOG_DBETA log(x)*(a - 1) + log(1 - x)*(b - 1) - lbeta(a, b)
 
-    if(give_log) {
-	return log(x)*(a - 1) + log(1 - x)*(b - 1) - lbeta(a, b);
-    }
+    if(give_log)
+	return R_LOG_DBETA;
+    else if (a + b >= xmax) /* beta(a,b) might be = 0 numerically */
+	return exp(R_LOG_DBETA);
     else {
 	double y;
 	y = beta(a, b);


------
Martin Maechler <maechler@stat.math.ethz.ch>
http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO D10	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><

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