beta_inc_inv sometimes not sufficiently accurate

[andreasnoack]

E.g.

julia> x
0.2142857142857142

julia> y = StatsFuns.beta_inc(1.5, 5.0, x)[1];

julia> (StatsFuns.beta_inc_inv(1.5, 5.0, y)[1] - x)/eps(x)
-140.0

I noticed this because R currently computes this more accurately.

I tried to look a bit into this and I noticed

SpecialFunctions.jl/src/beta_inc.jl

Line 974 in f2ccad8

iex = max(-5.0/a^2 - 1.0/p^0.2 - 13.0, -30.0)

. The expression in the first argument is not present in the publication that the algorithm is based and it's not in Applied Statistics Algorithms by Griffiths, P. and Hill, I.D but I do see the expression in http://lib.stat.cmu.edu/apstat/109 but without a comment about where it might come from. If I just set iex=-30 following the publications then I get

julia> (StatsFuns.beta_inc_inv(1.5, 5.0, y)[1] - x)/eps(x)
0.0

@simonbyrne Do you have any guess to where the extra expression that reduces iex come from?

[simonbyrne]

Ah, that's interesting, I hadn't notice the discrepancy from the paper. I have no idea!

[simonbyrne]

This seems to be the most recent work in the area https://doi.org/10.1007/s11075-016-0139-2 / https://arxiv.org/abs/1605.03503v1

[simonbyrne]

A relatively recent algorithm was given in [5].

[5] is from 1992.

[andreasnoack]

Thanks for the link. It looks like something that would interesting to try out instead of the current implementation. The paper also solved the mystery here since it cites the remark that introduces the adjustment of ACU. It's https://rss.onlinelibrary.wiley.com/doi/epdf/10.2307/2347779. However, they provide no motivation for the coefficients. The original implementations from Applied Statistics are all single precision so maybe the formula requires some adjustment for double precision but that it's impossible to tell from the comment.