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 …

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.473779306…

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…

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 t…

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::Bi…

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 re…

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

miguelraz April 28, 2021, 1:05am 2

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

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

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

Related topics