Codecov Report

Attention: Patch coverage is 80.09479% with 84 lines in your changes missing coverage. Please review.

Project coverage is 85.10%. Comparing base (3de6038) to head (bbf5468).

@@            Coverage Diff             @@
##           master    #1905      +/-   ##
==========================================
- Coverage   85.99%   85.10%   -0.90%     
==========================================
  Files         144      144              
  Lines        8685     8774      +89     
==========================================
- Hits         7469     7467       -2     
- Misses       1216     1307      +91     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

Given that, according to the docstring

#1907 is related.

#1907 is related.

17154a2 fixes it.

My personal opinion:
rand(d::Distribution) should remain a valid method (I don't think anyone is trying to change this), and its return type should be predictable by knowing what d is. It would also be useful to have a method that tells you what type rand(d::Distribution) returns without running rand(d::Distribution). eltype(d::Distribution) is logical for this.

This is orthogonal to the fact that Distributions should probably support rand(d::Distribution, T::Type).

In an ideal situation, I think we would define rand(d::Distribution, T::Type) for distributions (to allow for potential type-specific optimizations and evade unnecessary typecasts) and then dispatch rand(d::Distribution) = rand(d, eltype(d)).

its return type should be predictable by knowing what d is

I became convinced that in general this is not possible. It basically means trying to manually figure out type inference, a task that even the Julia compiler can fail at. Of course, for some number types and distributions it can be done (as also the compiler will succeed in many instances). But in general it's far from trivial. As in the initial stages of Distributions, one could be restrictive and limit samples to Int and Float but that's IMO not desirable either. Therefore I think it's best to remove eltype (which IMO has also been a weird name since I don't view distributions as containers or collections) from the API. IMO it's just too brittle and currently to inconsistent. Of course, that doesn't prevent users from trying to figure out the (element) type of samples with eg the help of the compiler, it's just not an official feature anymore.

A few distributions still give me Float64 after this PR, while others work fine:

julia> dist = Semicircle(50.0f0)
Semicircle{Float32}(r=50.0f0)

julia> rand(dist, 5)
5-element Array{Float64, 1}:
  36.80671953487414
 -18.355635129900335
 -11.701855436648922
 -21.444118928985656
  -5.80120463505302

julia> dist = JohnsonSU(0.0f0, 1.0f0, 0.0f0, 1.0f0)
JohnsonSU{Float32}(ξ=0.0f0, λ=1.0f0, γ=0.0f0, δ=1.0f0)

julia> rand(dist, 5)
5-element Array{Float64, 1}:
  0.5311696707298299
  1.632313034117999
  0.04951771555318912
  0.4721610259428258
 -3.052321854866766

julia> dist = Chisq(5.0f0)
Chisq{Float32}(ν=5.0f0)

julia> rand(dist, 5)
5-element Array{Float32, 1}:
 15.465032
 1.888659
 7.013455
 4.258529
 3.9611576

Can this be fixed?

The reason is that for quite a few of the older, probably less used and definitely less frequently updated distributions the rand implementation contains hardcoded Float64, either as literals or as calls of rand(rng) or randn(rng). I fixed the SemicircleandJohnsonSU`, but I'm not sure if it's possible (or even desirable) to address all of these in a single PR (it's already a quite massive undertaking). Maybe it would be best to open issues for others and address them in follow-up PRs (most (all?) other such changes in this PR are included to fix open issues).

This is not the first time I have went on a rant about how Distributions in this package are effectively unordered containers with infinite size, but I have stumbled upon new ammunition:
https://docs.julialang.org/en/v1/stdlib/Random/#A-simple-sampler-without-pre-computed-data

This section of the Julia stdlib documentation claims:

Given a collection type S, it's currently assumed that if rand(::S) is defined, an object of type eltype(S) will be produced. In the last example, a Vector{Any} is produced; the reason is that eltype(Die) == Any. The remedy is to define Base.eltype(::Type{Die}) = Int.

This led me to JuliaLang/julia#31787 and JuliaLang/julia#27756

I believe that it is actually deeper Random.jl precedent that typeof(rand(s)) == eltype(s). Do we really want to break this precedent in Distributions.jl?

Additional rambling

On a similar note, I do not think that it is controversial that eltype(d) should also be the type returned by minimum(d), maximum(d), mean(d), mode(d), etc., following Base precedent.

I also believe that it is reasonable that quantile(d, x) should return eltype(d); it is analogous to getindex to some extent. When would it make sense for rand and quantile to return different types? In other words, I believe that the following behavior is not ideal:

julia> d = Normal(0f0, 1f0)
Normal{Float32}(μ=0.0f0, σ=1.0f0)

julia> typeof(quantile(d, .5)) == typeof(rand(d))
false

For many distribution types, perhaps the real conversation we should be having is if we should force partype(D == eltype(D). For distributions parameterized by their extrema, e.g. uniform, I think this is very reasonable. Uniform(0f0, 1f0) is basically saying "construct an infinite unordered container with a uniform distribution with minimum == 0f0, maximum == 1f0."

For distributions parameterized by mean/std, similar logic definitely holds. It is not a stretch to extend this for other location/spread parameters. Shape parameters are ambiguous, but if someone types Beta(1f0, 1f0), I would suspect that they have a good reason for it.

For distributions like Binomial, I would argue that the type of n is more useful than the type of p. Unfortunately, the type of n is currently hardcoded as Int while the type of p is user-defined.

Yet further rambling

cdf is a strange cousin of indexin and findfirst, which I believe are guaranteed to return keytype. Actually implementing an explicit keytype(d) == Float64) may be a reach though, as I suppose this also implies getindex support, which quantile is not quite identical to.