Hi everyone, I was comparing a colleague’s MATLAB code to my version (we are coding in parallel, and I am using Julia. I am also trying to convince him to move to Julia. Making progress… :slight_smile: ), and I noticed that the beta function returns something different for large x or y. Because of some other related functions these values are fed into, there are some differences that emerge in our results. I elaborate on this in the following code: using SpecialFunctions beta1 = 3.6e13; alpha1 = 0.563483398; a13 = -1.27e29; julia> beta1 3.6e13 julia> alpha1 0.563483398 julia> beta(beta1,1/alpha1+1) 3.9891491767885927e-38 # MATLAB's version >> beta(beta1,1/alpha1+1) ans = 4.47377930618112e-38 The differences emerge in our codes because we use the following expression, which result in different answers. julia> a13 -1.27e29 julia> a13*beta1*beta(2*beta1,1/alpha1+1) -26651.76910015019 # MATLAB's version >> a13*beta1*beta(2*beta1,1/alpha1+1) ans = -21558.4828041656 Would anyone have any insight about why these two values differ? I tried reading up on this and found some useful issues like #14349 and #14620 but I am wondering if someone has any additional comments or insights regarding what is likely happening here.

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

It’s impossible to help without seeing your full code: Please read: make it easier to help you

The full code is too long to post. However, I think what I have posted highlights the difference and provides all the values that are needed to compute the expressions. If you copy-paste what I have posted, it should run on your machine. The beta function is in SpecialFunctions.jl I believe.

Ah, I see the SpecialFunctions.jl import is now at the top of the file. That should do it. I don’t know enough about the floating point specifics of the beta function to tell you the correct answer, but you can always use BigFloats to check your answer in Julia: julia> a13*beta1*beta(big(2*beta1),big(1)/alpha1+1) -26651.76910015004456135864648326051363960481520541247651058910044861872889956144

I think the question is more about why these values differ so much between MATLAB and Julia. I’m not an expert in floating point arithmetic, so maybe someone with more knowledge can help.

This is weird [image] miguelraz: julia> a13*beta1*beta(big(2*beta1),big(1)/alpha1+1) -26651.76910015004456135864648326051363960481520541247651058910044861872889956144 I get NaN here. Julia 1.6.0, SpecialFunctions 1.3.0

I think Julia’s answer is closer to correct, but with this large of a number it’s easy to blow up or hit a catastrophic cancellation. To take this back one step, a naive definition for beta is \beta(p, q) = \frac{\Gamma(p) \Gamma(q)}{\Gamma(p+q)} You wouldn’t want to use this implementation in the real world, though, because those intermediate values are huge. Like 10^{472392288679098} huge. So I’m quite certain this isn’t the implementation in either language. julia> gamma(3.6e13) Inf But, hey, we can opt into higher precision floating point here: julia> gamma(big(3.6e13)) 7.882861104028365451673058070377064242970022825665897715265241972922672893761775e+472392288679098 julia> p = big(3.6e13) 3.6e+13 julia> q = 1/big(0.563483398)+1 2.774675178628776489804929213602311938188089959295862682297487944939957397047059 julia> gamma(p)*gamma(q)/gamma(p+q) 3.989149176788555366691087811957298867890005562542321569109120632926257477147039e-38 julia> beta(p, q) 3.989149176788555366691087811957298867890005562542321569109120643773484257104723e-38 So I’m not a numericist, but I think this is looking pretty good.

This is very helpful! And, another reason to hopefully convince him to switch to Julia soon :wink:

[image]gamma(::BigFloat) overflows on Windows @which gamma(big(3.6e13)) Weird. From my side julia> using SpecialFunctions julia> gamma(big(3.6e13)) 7.882861104028365451673058070377064242970022825665897715265241972922672893761775e+472392288679098 julia> @which gamma(big(3.6e13)) gamma(x::BigFloat) in SpecialFunctions at /Users/gerzaguet/.julia/packages/SpecialFunctions/mFAQ4/src/gamma.jl:578 (testSpecialFunctions) pkg> st Project testSpecialFunctions v0.1.0 Status `~/testSpecialFunctions/Project.toml` [276daf66] SpecialFunctions v1.3.0 (testSpecialFunctions) pkg>

Mathematica also likes Julia’s answer: In[9]:= x=36*10^12; y=1+10^9/563483398; Beta[N[x,30],N[y,30]] Out[10]= 3.98914917678855*10^-38 The fact that the number of accurate digits (in MMA’s estimation) drops from 30 to about 15 suggests a condition number of about 10^15, which is very nearly the reciprocal of double precision roundoff. In fact, if I tell MMA to use 16 digits in the parameters, it reports the result as simply 4*10^(-38). Thus both Julia and MATLAB are reporting the same, right answer, but you can’t trust more than one significant digit due to the limit of double precision. I suspect that Julia’s library may use extended precision to get more accurate digits than you have a right to expect. That is, it’s treating your parameters as being given with a lot more precision than they actually have, as BigFloats do. EDIT: In fact, if beta1 and/or alpha1 has actually been rounded off to 9 significant digits as it would appear, then the result is utterly meaningless. The next, unknown digit could change the answer by more than 100%.

Mathematica seems to agree with MATLAB here: [image]

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

General Usage

Topic		Replies	Views
gamma(::BigFloat) overflows on Windows General Usage specialfunctions	7	1184	April 29, 2021
Result of exponentials are (sometimes) different from Matlab Numerics	5	529	February 15, 2024
Wrong results in calculation Performance	21	3601	March 12, 2019
Beta Distribution - full precision may not have been achieved Statistics	7	2038	June 13, 2018
Unexpected result about integer General Usage	13	1623	August 17, 2018

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

Related topics