Integrating beta-function-like function over compact support [0,alpha], alpha < 1

Question

I would like to integrate a following function named betalog

 g <- function(x,a,b){
         if (a < 0 | b < 0) stop()
         temp <- (a-1)*log(x) + (b-1)*log(1-x)
         return( exp(temp) )
       }

 betalog<- function(x,a,b)
   {
    temp <- g(x=x,a=a,b=b)* log(x/(1-x))
    return( temp )
   }

The function g is integrand of the beta function. In theory, betalog should be integrable over any [0,alpha] interval if 0 < alpha < 1, and a > 0, b >0. However, I cannot numerically integrate betalog with very small a:

 a <- 0.00001
 b <- 1
 alpha <- 0.5
 integrate(betalog,a=a,b=b,lower=0,upper=alpha,subdivisions=1000000L)

 Error in integrate(betalog, a = a, b = b, lower = 0, upper = alpha, subdivisions =  
 1000000L) : 
 non-finite function value

In fact, I cannot even compute the incomplete beta function using R integrate function when a is very small:

 integrate(g,a=a,b=b,lower=0,upper=alpha,subdivisions=1000000L)

 Error in integrate(g, a = a, b = b, lower = 0, upper = alpha, subdivisions = 1000000L) : 
 roundoff error is detected in the extrapolation table

Can anyone gives me tip to integrate such incomplete beta-like function in R?

Answer 1

> betalog(0, a, b)
[1] -Inf

Your function is singular at the lower bound. Recall that to compute an improper integral you must replace the singular bounds with dummy variables and take the limit from the correct side towards that bound. In particular,

 >  integrate(betalog,a=a,b=b,lower=0.000001,upper=alpha,subdivisions=10000000L)
 -94.60292 with absolute error < 0.00014
 >  integrate(betalog,a=a,b=b,lower=.Machine$double.xmin * 1000,upper=alpha,subdivisions=10000
 -244894.7 with absolute error < 10
 > integrate(betalog,a=a,b=b,lower=.Machine$double.xmin,upper=alpha,subdivisions=10000000L)
 Error in integrate(betalog, a = a, b = b, lower = .Machine$double.xmin,  :
   non-finite function value

I suspect that your integral diverges, but this might be tricky since even state-of-the-art symbolic algebra systems can't prove that:

http://www.wolframalpha.com/input/?i=Integral%28x%5E%280.00001+-1%29+ln%28x%2F%281-x%29%29%2C+x%2C0%2C+0.5%29

Whatever the case, R is not the correct tool for this problem.