# Accurate implementation of the logarithm of the incomplete Beta function in C++?

I need an accurate implementation (for use in C/C++) of the logarithm of the incomplete Beta function:

$$\log \mathrm{B}(x,y;\alpha,\beta) = \log \int_x^y t^{\alpha-1}(1-t)^{\beta-1}\mathrm{d}t$$

with domain $\alpha,\beta>0$, $0\le x\le y \le 1$. I need this to avoid overflows with very large values of the incomplete Beta function.

Boost and GSL provide direct implementations of $\mathrm{B}(x,y;\alpha,\beta)$ or its regularized version. This won't work for me, because from these there's no way to get the value of $\log \mathrm{B}(x,y;\alpha,\beta)$ without going through the large values I am trying to avoid (or at least I don't know how).

Update: I found that R has an implementation of the logarithm of the incomplete Beta function (see https://stat.ethz.ch/R-manual/R-patched/library/stats/html/Beta.html). Can I use this from C++? Will it be as fast as native C code?
• Have you seen this question/answer for the "lower" incomplete beta function, $B(0,y;a,b)$? – horchler May 23 '14 at 22:38
• Have you tried $\ln B(x,y;\alpha,\beta) = \ln B(0,y;\alpha,\beta) - \ln B(0,x;\alpha,\beta) = \ln \frac{B(0,y;\alpha,\beta)}{B(0,x;\alpha,\beta)}$? Or are the numerics not good (e.g., catastrophic cancellation)? – horchler Jun 5 '14 at 0:57
• @horchler Your equation is not correct. The correct equation is $\mathrm{B}(x,y;\alpha,\beta) = \mathrm{B}(y;\alpha,\beta) - \mathrm{B}(x;\alpha,\beta)$, without logarithms. – becko Jun 5 '14 at 12:31
Boost provides all the functions you need, so it should be sufficient. Instead of using the direct integral representation of the incomplete beta function, it's better to use the regularized incomplete beta function (see DLMF 8.17): $$B(x,y;a,b) = (I_y(a,b)-I_x(a,b))B(a,b).$$ This way, $I_x(a,b)=B(x;a,b)/B(a,b)$ will have values in $[0,1]$, so that $$\log B(x,y;a,b) = \log(I_y(a,b)-I_x(a,b)) + \log B(a,b).$$ is much easier to compute. This way the only source of numerical instability, assuming the underlying functions are stable, would be the subtraction when $x$ and $y$ are very close; in that case, if it ever arises, the direct integral representation might work.
The regularized incomplete beta function is given by ibeta in Boost, and the log-Beta function can be calculated using the log-gamma function lgamma, also in Boost. Boost's implementation of $I_x$ is based on the paper Algorithm 708: Significant digit computation of the incomplete beta function ratios by Didonato and Morris.