Cannot sample from legal NegativeBinomial distribution

[damonbayer]

I have noticed that some arguments allowed in the NegativeBinomial distribution lead to distributions that cannot be sampled from. Here are two examples:

julia> r_1 = nextfloat(0.0)
5.0e-324

julia> p_1 = 1.0
1.0

julia> allowed_nb_1 = NegativeBinomial(r_1, p_1)
NegativeBinomial{Float64}(r=5.0e-324, p=1.0)

julia> rand(allowed_nb_1)
ERROR: LoadError: DomainError with 0.0:
Gamma: the condition θ > zero(θ) is not satisfied.

# This happens because sampling is done by rand(rng, Poisson(rand(rng, Gamma(d.r, (1 - d.p)/d.p))))

julia> Gamma(r_1, (1 - p_1)/p_1)
ERROR: LoadError: DomainError with 0.0:
Gamma: the condition θ > zero(θ) is not satisfied.

julia> r_2 = nextfloat(0.0)
5.0e-324

julia> p_2 = nextfloat(0.0)
5.0e-324

julia> allowed_nb_2 = NegativeBinomial(r_2, p_2)
NegativeBinomial{Float64}(r=5.0e-324, p=5.0e-324)

julia> rand(allowed_nb_2)
ERROR: LoadError: DomainError with NaN:
Poisson: the condition λ >= zero(λ) is not satisfied.

# This happens because sampling is done by rand(rng, Poisson(rand(rng, Gamma(d.r, (1 - d.p)/d.p))))

julia> allowed_g_2 = Gamma(r_2, (1 - p_2)/p_2)
Gamma{Float64}(α=5.0e-324, θ=Inf)

julia> rand(allowed_g_2)
NaN

[simsurace]

The first case can certainly be fixed by a separate branch in rand(::NegativeBinomial) if isone(p), but I'm not sure what should be done in the second one. What would you expect rand(Gamma(nextfloat(0.0), Inf)) to return? rand(Poisson(Inf)) will throw anyways, and we don't want to return Inf because that will make the function type-unstable. Also, typemax(Int) is much smaller than prevfloat(Inf).

[devmotion]

I think the second case, and generally such issues with unbounded discrete distributions where Int is not sufficient, would be fixed by https://github.com/JuliaStats/Distributions.jl/pull/1433.