3
$\begingroup$

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).

Any advice?

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?

|cite|improve this question
$\endgroup$
  • 1
    $\begingroup$ Have you seen this question/answer for the "lower" incomplete beta function, $B(0,y;a,b)$? $\endgroup$ – horchler May 23 '14 at 22:38
  • $\begingroup$ +1 @horchler Thanks. Any idea of how I can use that to get the logarithm of general incomplete beta function? $\endgroup$ – becko Jun 4 '14 at 12:24
  • $\begingroup$ 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)? $\endgroup$ – horchler Jun 5 '14 at 0:57
  • $\begingroup$ @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. $\endgroup$ – becko Jun 5 '14 at 12:31
3
$\begingroup$

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.

|cite|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.