Differences in MATLAB's versus Julia's beta(x, y) for large x or y?

General Usage
18

To build on other’s answers that express beta() as \Gamma(p)\Gamma(q)/\Gamma(p+q), it seems you can trivially reproduce Matlab’s behavior in Julia using loggamma():

julia> p = beta1
3.6e13

julia> q = 1/alpha1 + 1
2.7746751786287764

julia> exp(loggamma(p) + loggamma(q) - loggamma(p+q))
4.4737793061811207e-38

julia> exp(loggamma(big(p)) + loggamma(big(q)) - loggamma(big(p)+q))
3.98914917678856144745513696868341367384631253052121882201338979974154777667653e-38

I’d bet on round-off error, as others have pointed out.

EDIT: loggamma() is preferred over lgamma().

6 Likes