Logarithm of incomplete Beta function for large $\alpha,\beta$ R's function pbeta is supposed to calculate the incomplete regularized Beta function. If the flag log=TRUE is passed as an argument, the logarithm of the result is returned instead.
I found the following numerical error. Evaluating:
pbeta(0.5555555, 1925.74, 33.7179, log=TRUE)

returns -Inf, which is obviously wrong. The correct result is (calculated using Mathematica):
$$\log\mathrm{I}_{0.5555555}(1925.74,33.7179) = -994.767594138466967$$
For all the other values I have tested, pbeta gives the correct answer. I've inspected the source code of pbeta but cannot pin-point the error. 
I need an accurate method to compute the logarithm of the incomplete Beta function for large values of $\alpha,\beta$ (on the order of thousands). Can someone suggest a fix on pbeta or a workaround? Or perhaps a different library?
R bug report: https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=16332
 A: The Beta incomplete function can be written with the help of the Gauss hypergeometric function:
$$B_x(a,b)=\frac{1}{a}x^a{(1-x)}^b F(a+b,1,a+1,x)$$
and the CDF of the $Beta(a,b)$ distribution evaluated at $x$ is 
$$B_x(a,b)/B(a,b).$$
A good implementation of the Gauss hypergeometric function is provided in the gsl package - a wrapper for the special functions of the Gnu Scientific Library. So you can write the log-CDF like this:
library(gsl)
logpbeta <- function(x,a,b) log(hyperg_2F1(a+b,1,a+1,x)) + a*log(x)+b*log(1-x)-log(a) - lbeta(a,b)

And it gives the same result as Mathematica for your example:
> logpbeta(0.5555555, 1925.74, 33.7179)
[1] -994.7676

I don't know for which range the parameters the result is correct. And I have not find the answer in the GNU reference manual. 
Note that the zero you obtain with pbeta seems to be due to the evaluation of $B_x(a,b)$ followed by the log-transformation:
> x <- 0.5555555 
> a <- 1925.74
> b <- 33.7179
> log(hyperg_2F1(a+b,1,a+1,x)*x^a*(1-x)^b/a)
[1] -Inf
> hyperg_2F1(a+b,1,a+1,x)
[1] 2.298761
> x^a*(1-x)^b/a
[1] 0

