R devel - Sep 2013 - dbeta may hang R session for very large values of the shape parameters

Dear all,

we received a bug report for betareg, that in some cases the optim call in
betareg.fit would hang the R session and the command cannot be interrupted by
Ctrl-C?

We narrowed down the problem to the dbeta function which is used for the log
likelihood evaluation in betareg.fit.

Particularly, the following command hangs the R session to a 100% CPU usage in
all systems we tried it (OS X 10.8.4, Debian GNU Linux, Ubuntu 12.04) with
either R-3.0.1 and with the R-devel version (in all systems I waited 3 minutes
before I kill R):

## Warning: this will hang the R session
dbeta(0.9, 1e+308, 10)

Furthermore, through a trial and error investigation using the following code

## Warning: this will hang the R session
x <- 0.9
for (i in 0:100) {
 a <- 1e+280*2^i
 b <- 10
 cat("shape1 =", a, "\n")
 cat("shape2 =", b, "\n")
 cat("Beta density", dbeta(x, shape1 = a, shape2 = b), "\n")
 cat("===\n")
}

I noticed that:
* this seems to happen when shape1 is about 1e+308, seemingly irrespective of
the value of shape2 (run the above with another value of b), and as it appears
only when x>=0.9 and x < 1 (run the above lines with x <- 0.89999 for
example and everything works as expected).
* similar problems are encountered for small x values when shape2 is massive.

I am not sure why this happens but it looks deep to me. The temporary fix for
the purposes of betareg was a hack (a simple if command that returns NA for the
log likelihood if any shape parameter has values greater than 1e+300 say).

Nevertheless, I thought that this is an issue worth reporting to R-devel
(instead of R-help), especially since dbeta may be used within generic
optimisers and figuring that dbeta is the problem can be hard --- it took us
some time before we started suspecting dbeta.

Interestingly, this appears to happen close to what R considers infinity. Typing
1.799e+308
into R returns Inf.

I hope the above limited in scope analysis is informative.

Best regards,
Ioannis



-- 
Dr Ioannis Kosmidis
Department of Statistical  Science,
University College,
London, WC1E 6BT, UK
Webpage: http://www.ucl.ac.uk/~ucakiko

The issue is underflow to zero in bd0 (C file src/nmath/bd.c).  We'll 
fix that, but given the R 3.0.2 RC is in code freeze and this has 
existed for years, not for 3.0.2.

On 18/09/2013 23:52, Kosmidis, Ioannis wrote:
> Dear all,
>
> we received a bug report for betareg, that in some cases the optim call in
betareg.fit would hang the R session and the command cannot be interrupted by
Ctrl-C?
>
> We narrowed down the problem to the dbeta function which is used for the
log likelihood evaluation in betareg.fit.
>
> Particularly, the following command hangs the R session to a 100% CPU usage
in all systems we tried it (OS X 10.8.4, Debian GNU Linux, Ubuntu 12.04) with
either R-3.0.1 and with the R-devel version (in all systems I waited 3 minutes
before I kill R):
>
> ## Warning: this will hang the R session
> dbeta(0.9, 1e+308, 10)
>
> Furthermore, through a trial and error investigation using the following
code
>
> ## Warning: this will hang the R session
> x <- 0.9
> for (i in 0:100) {
>   a <- 1e+280*2^i
>   b <- 10
>   cat("shape1 =", a, "\n")
>   cat("shape2 =", b, "\n")
>   cat("Beta density", dbeta(x, shape1 = a, shape2 = b),
"\n")
>   cat("===\n")
> }
>
> I noticed that:
> * this seems to happen when shape1 is about 1e+308, seemingly irrespective
of the value of shape2 (run the above with another value of b), and as it
appears only when x>=0.9 and x < 1 (run the above lines with x <-
0.89999 for example and everything works as expected).
> * similar problems are encountered for small x values when shape2 is
massive.
>
> I am not sure why this happens but it looks deep to me. The temporary fix
for the purposes of betareg was a hack (a simple if command that returns NA for
the log likelihood if any shape parameter has values greater than 1e+300 say).
>
> Nevertheless, I thought that this is an issue worth reporting to R-devel
(instead of R-help), especially since dbeta may be used within generic
optimisers and figuring that dbeta is the problem can be hard --- it took us
some time before we started suspecting dbeta.
>
> Interestingly, this appears to happen close to what R considers infinity.
Typing
> 1.799e+308
> into R returns Inf.
>
> I hope the above limited in scope analysis is informative.
>
> Best regards,
> Ioannis
>
>
>

-- 
Brian D. Ripley,                  ripley at 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 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595