Distributions uses an old version of Rmath

[StatsMan1953]

R v-4.1.3 (released 10 March 2022) improved accuracy of distribution quantile functions, but Distributions continues to use an older version of Rmath.

Example:
quantile(Beta(0.01,0.1),0.92)
returns 1.0 instead of 0.780301421091985
Older versions of R at least threw a warning "full precision may not have been achieved in 'qbeta'", but Julia gives no indication that the result is completely wrong.

[devmotion]

Thanks for opening the issue!

This problem is interesting in so far as it concerns two different packages but unfortunately it seems we can't address it in Distributions. Some explanations:

• Distributions just uses the implementations in StatsFuns for evaluating certain functions such as the quantiles of Beta distributions (and some other distributions as well), i.e., it just defines quantile(d::Beta, p::Real) = StatsFuns.betainvcdf(d.alpha, d.beta, p)
• StatsFuns.jl < 1 used Rmath.jl for many evaluations, including betainvcdf:

  julia> betainvcdf(0.01, 0.1, 0.92)
  full precision may not have been achieved in 'qbeta'
  1.0

• StatsFuns.jl >= 1.0 uses native Julia implementations based on SpecialFunctions.jl and HypergeometricFunctions.jl for many functions, including betainvcdf:

  julia> betainvcdf(0.01, 0.1, 0.92)
  1.0

• Recent versions of Rmath.jl use Rmath_jll which is built from https://github.com/JuliaStats/Rmath-julia
• Rmath-julia is outdated it seems (but in StatsFuns >= 1 this is not relevant for the problem here): see e.g. Standalone Rmath Version bump: 4.0.4 -> 4.1.0 Rmath-julia#42

I'll transfer the issue to StatsFuns.

[devmotion]

Maybe the issue should actually be transferred to SpecialFunctions since betainvcdf just calls SpecialFunctions.beta_inc_inv:

StatsFuns.jl/src/distrs/beta.jl

Line 83 in e5699db

betainvcdf(α::Real, β::Real, p::Real) = first(beta_inc_inv(α, β, p))

[devmotion]

Indeed, it's a SpecialFunctions issue. With SpecialFunctions 2.1.4 I get

julia> using SpecialFunctions

julia> beta_inc_inv(0.01, 0.1, 0.92)
(0.7803014210919854, 0.21969857890801453)

Seems one of the fixes in JuliaMath/SpecialFunctions.jl#399 accidentally broke this example here.

[devmotion]

With JuliaMath/SpecialFunctions.jl#401 I get

julia> betainvcdf(0.01, 0.1, 0.92)
0.7803014210919854

as before with SpecialFunctions 2.1.4. This seems reasonably accurate:

➜ ~ wolframscript -c "N[InverseCDF[BetaDistribution[1/100, 1/10], 92/100]]"
0.7803014210919872

[devmotion]

I'll close this issue, the problem is fixed in SpecialFunctions 2.1.6. Upgrading your packages should resolve the issue.