[Rd] dbeta may hang R session for very large values of the shape parameters

[Ioannis Kosmidis]

Dear Brian,

many thanks, that is great.

Best regards,
Ioannis

On 19/09/13 12:15, Prof Brian Ripley wrote:
> 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
>>
>>
>>
>
>



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