fix type stability of sampling from Chisq, TDist, Gamma

[Red-Portal]

This addresses #1884

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 85.99%. Comparing base (13029c0) to head (2091cd7).
Report is 2 commits behind head on master.

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1885   +/-   ##
=======================================
  Coverage   85.99%   85.99%           
=======================================
  Files         144      144           
  Lines        8666     8666           
=======================================
  Hits         7452     7452           
  Misses       1214     1214           

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

[quildtide]

Should probably add method eltype(d::Gamma) = partype(d) and friends.

There are many side effects associated with changing Gamma. The following list is not exhaustive, but I think they will be impacted:

• Beta
• BetaPrime
• InverseGamma
• Chi

Chisq will have some side effects too, e.g.:

• FDist

[Red-Portal]

Hi @quildtide ,

Okay so to be precise, I'll do the following:

• Identify all the distributions that are affected by sampling Gammas. If the type of the sampled gamma affects the type of the final sample, add eltype = partype
• Add tests for the type stability

Would that be sufficient?

[quildtide]

Not a maintainer, so don't take my word for granted. If the route this pull request takes is chosen, then the things you propose doing would indeed be a good idea.

Beta has two sampling paths, IIRC, where it sometimes uses a Gamma sampler, and sometimes uses a different sampler, so you will probably need to alter the other path too, if the Gamma path is returning partype values now (I think it is).

---

I know that some of the maintainers have pushed back to some degree in the past on making eltype == partype for Distributions that aren't already doing that (Normal and a few friends), but @devmotion seems to be calling the shots right now, and he seems to back the core idea of this pull req for TDist.

[devmotion]

This is a slightly different case and easily broken (e.g., when T is not a floating point number type). In the TDist and Gamma case we just try to avoid promotions of a sample from another rand call, whereas this case goes deeper into the question of how rand should behave wrt parameters etc. (see also #1433 (comment)).

Suggested change

	rand(rng::AbstractRNG, d::Exponential{T}) where {T} = xval(d, randexp(rng, T))
	rand(rng::AbstractRNG, d::Exponential) = xval(d, randexp(rng))

[Red-Portal]

If we do this though, the return type of rand(Gamma(Float32, Float32)) changes depending on the value of the shape parameter because shape == 1 samples from Exponential. (This is why the tests are currently failing.) Should we let this happen? I imagine some people will be super surprised by such behavior.

[Red-Portal]

Hi @devmotion could you comment on this?

[devmotion]

That's a good point. I didn't realize that GammaMTSampler already respects the parameter types (but samples are not necessarily of the parameter type:

Distributions.jl/src/samplers/gamma.jl

Lines 175 to 182 in 13029c0

	d = shape(g) - 1//3
	c = inv(3 * sqrt(d))
	# Pre-compute scaling factor
	κ = d * scale(g)
	# We also pre-compute the factor in the squeeze function
	return GammaMTSampler(promote(d, c, κ, 331//10_000)...)

).

[devmotion]

I think that's an argument for using the same approach as for Normal here, until we move to a better/different API:

Suggested change

	rand(rng::AbstractRNG, d::Exponential{T}) where {T} = xval(d, randexp(rng, T))
	rand(rng::AbstractRNG, d::Exponential{T}) where {T} = xval(d, randexp(rng, float(T)))

[devmotion]

Can you also add a test for Exponential?

[Red-Portal]

• I also added back the type stability tests for Chisq since we now know that the gamma samplers respect the parameter type.
• I removed the type stability test for Entropy(TDist) because it is not type stable (it returns Union{Float32, Float64}) because $\nu = \mathrm{Inf}$ returns entropy(Normal), which is also not type stable.