Regression v0.9.15 to v0.9.16: betainvcdf stalls for small values

[oschulz]

With StatsFuns v0.9.16 this doesn't terminate (at least not within a few minutes):

julia> StatsFuns.betainvcdf(1.0, 1.0, 1e-21)

StatsFuns v0.9.15 gives the correct result:

julia> StatsFuns.betainvcdf(1.0, 1.0, 1e-21)
1.0e-21

julia> versioninfo()
Julia Version 1.7.2
Commit bf53498635 (2022-02-06 15:21 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) Gold 6140 CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, skylake-avx512)

[devmotion]

Seems like another regression due to the switch from Rmath to SpecialFunctions for some functions. (Unfortunate for users but it's astonishing to see how many bugs for special cases are discovered in SpecialFunctions now that it is used more intensively by StatsFuns and Distributions.)

[oschulz]

I did take me a while to track it down (happened deep in a sampler in BAT.jl), but as far as I'm concerned getting rid of Rmath is so worth it. :-)

[oschulz]

Seems to happen only for d.α ≈ 1 && d.β ≈ 1 with u < 1e-19 or so.

[andreasnoack]

Do you have an example where one or both of the parameters are different from one. It looks like the algorithm has an issue specifically when they are one. I think we should probably just handle those as special cases as the quantile function also has a closed form in those cases.

[oschulz]

Do you have an example where one or both of the parameters are different from one.

I tried a few things, but couldn't find any other case that stalls - the problem only seems to occur if both α ≈ 1 and β ≈ 1.

[andreasnoack]

You are using approx here again. I'd be interested an any non-unit values of a and b that trigger the issue. As I read the code, it's only exact unity that causes an issue so it would be good to know if there are other cases as well.

[oschulz]

Ah, sorry, I misunderstood. Yes, there seem to be cases that require approx (maybe only for the second argument):

StatsFuns.betainvcdf(1.0, 1.000000000001, 1e-21)

also stalls, though

StatsFuns.betainvcdf(1.000000000001, 1.0, 1e-21)

does not.

[andreasnoack]

This will be fixed by https://github.com/JuliaMath/SpecialFunctions.jl/pull/393

[oschulz]

Thanks!!

[andreasnoack]

Should be fixed by https://github.com/JuliaRegistries/General/pull/55444. Please report back if you experience more regressions.

[oschulz]

@andreasnoack , @devmotion can we reopen this? I ran into another case that stalls:

betainvcdf(2.0, 1.0, 2.0e-280)

Interestingly, betainvcdf(2.0, 1.0, 2.0e-276) and betainvcdf(2.0, 1.0, 2.0e-281) do not stall.

[oschulz]

Thanks @devmotion .

[andreasnoack]

So I just tried the Fortran version from statlib and it has the same issue with these inputs so at least it isn't a bug that we introduced. I'll try to figure out if something can be improved wrt precision here to avoid the infinite loop or if we simply need to bound the number of iterations.

[oschulz]

StatsFuns v0.9.15 runs betainvcdf(2.0, 1.0, 2.0e-280) without problems, but that used a different lib, right?

[andreasnoack]

StatsFuns v0.9.15 runs betainvcdf(2.0, 1.0, 2.0e-280) without problems, but that used a different lib, right?

Yeah. It's RMath based. It probably originates from the same Fortran source but they have then most likely applied modifications to the ones we have to go through here.

https://github.com/JuliaMath/SpecialFunctions.jl/pull/396 should fix the immediate issue here.

[oschulz]

Thanks!

[devmotion]

All examples mentioned in this issue work without hang with SpecialFunctions 2.1.5:

julia> betainvcdf(1.0, 1.0, 1e-21)
9.99999999999998e-22

julia> betainvcdf(1.0, 1.000000000001, 1e-21)
9.999999999989961e-22

julia> betainvcdf(1.000000000001, 1.0, 1e-21)
1.0000000000483575e-21

julia> betainvcdf(2.0, 1.0, 2.0e-276)
1.4142135623731116e-138

julia> betainvcdf(2.0, 1.0, 2.0e-280)
1.4142135623731337e-140

julia> betainvcdf(2.0, 1.0, 2.0e-281)
4.4721359549997026e-141

Out of curiosity I also compared the results with Mathematica:

➜ ~ wolframscript -c "InverseCDF[BetaDistribution[1, 1], 10^-21] // N"
1.*^-21
➜ ~ wolframscript -c "InverseCDF[BetaDistribution[1, 1000000000001/1000000000000], 10^-21] // N"
9.999999999989997*^-22
➜ ~ wolframscript -c "InverseCDF[BetaDistribution[1000000000001/1000000000000, 1], 10^-21] // N"
1.0000000000483612*^-21
➜ ~ wolframscript -c "InverseCDF[BetaDistribution[2, 1], 2*10^-276] // N"
1.414213562373095*^-138
➜ ~ wolframscript -c "InverseCDF[BetaDistribution[2, 1], 2*10^-280] // N"
1.4142135623730948*^-140
➜ ~ wolframscript -c "InverseCDF[BetaDistribution[2, 1], 2*10^-281] // N"
4.472135954999579*^-141

[oschulz]

Thanks @devmotion !