Sampling from `StudentizedRange` returns `NaN` for `k < 2`

[ForceBru]

In Distributions v0.25.120 rand(StudentizedRange(10, k)) seems to always produce NaN for k < 2:

julia> rand(StudentizedRange(10, 1.1), 10)'
1×10 adjoint(::Vector{Float64}) with eltype Float64:
 NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN
julia> rand(StudentizedRange(10, 1.9), 10)'
1×10 adjoint(::Vector{Float64}) with eltype Float64:
 NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN
julia> rand(StudentizedRange(10, 2), 10)'
1×10 adjoint(::Vector{Float64}) with eltype Float64:
 0.237901  3.54309  0.883281  1.79658  0.70134  0.30904  0.692738  1.68206  1.65031  2.70354
julia> [mean(isnan, rand(StudentizedRange(10, k), 100)) for k in range(1.1, 2, 20)]'
1×20 adjoint(::Vector{Float64}) with eltype Float64:
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0

The last line shows that for 1.1 <= k < 2, the proportion of NaN in a sample is 1.0, so all samples are NaN.

[devmotion]

In Distributions v0.25.120

There's no recent changes to StudentizedRange in Distributions (last change of src/univariates/studentizedrange.jl was apparently 3 years ago), so I suspect the same issue has existed also in previous versions of Distributions?

Generally, samples from StudentizedRange are generated with inverse transform sampling using the qtukey quantile function from Rmath (exposed through StatsFuns). However, as reported in other places qtukey has difficulties for some parameter configurations and may return NaN instead of the correct values. Moreover, as noted in the R documentation,

A Legendre 16-point formula is used for the integral of ptukey. The computations are relatively expensive, especially for qtukey which uses a simple secant method for finding the inverse of ptukey. qtukey will be accurate to the 4th decimal place.

Maybe at this point we should just implement a pure Julia version of qtukey, drawing inspiration from SciPy.

---

Confirming that the parameter configurations above are problematic for qtukey:

julia> using Rmath

julia> all(isnan, qtukey.(0:0.01:1, 10, 1.1))
true

julia> all(isnan, qtukey.(0:0.01:1, 10, 1.9))
true

julia> !any(isnan, qtukey.(0:0.01:1, 10, 2))
true