Computes an approximation of the lbeta function which is equivalent to log(abs(Beta(a, b))) but avoids overflow by computing it with lgamma.
| 597 | // Computes an approximation of the lbeta function which is equivalent to |
| 598 | // log(abs(Beta(a, b))) but avoids overflow by computing it with lgamma. |
| 599 | static XlaOp Lbeta(XlaOp a, XlaOp b) { |
| 600 | // Beta(a, b) can be computed using Gamma as per |
| 601 | // http://dlmf.nist.gov/5.12.E1 as follows: |
| 602 | // Beta(a, b) = (Gamma(a) * Gamma(b)) / Gamma(a + b) |
| 603 | // |
| 604 | // To avoid overflow, we compute in the log domain. |
| 605 | // |
| 606 | // As per http://dlmf.nist.gov/4.8.E2 we can transform: |
| 607 | // Log(a * b) |
| 608 | // into: |
| 609 | // Log(a) + Log(b) |
| 610 | // |
| 611 | // Likewise, per https://dlmf.nist.gov/4.8.E4, we can turn: |
| 612 | // Log(a - b) |
| 613 | // into: |
| 614 | // Log(a) - Log(b) |
| 615 | // |
| 616 | // This means that we can compute Log(Beta(a, b)) by: |
| 617 | // Log(Gamma(a)) + Log(Gamma(b)) - Log(Gamma(a + b)) |
| 618 | return Lgamma(a) + Lgamma(b) - Lgamma(a + b); |
| 619 | } |
| 620 | |
| 621 | // Compute the Digamma function using Lanczos' approximation from "A Precision |
| 622 | // Approximation of the Gamma Function". SIAM Journal on Numerical Analysis |
no test coverage detected