Distribution sampling does not respect chosen number types

[ChrisRackauckas]

MWE:

data_dist = Beta(2.0f0, 2.0f0)
train_data = typeof(rand(data_dist, 1, 100)) # Matrix{Float64}

This is scary 😅

[devmotion]

There are multiple (open and closed) issues about this behaviour: e.g. https://github.com/JuliaStats/Distributions.jl/issues/1071 and linked issues.

Generally, the parameters are just treated as parameters, distinct from the type of the samples. The weird thing is Normal here which behaves differently (somehow it was changed at some point). I opened a PR a while ago that fixes these inconsistencies (https://github.com/JuliaStats/Distributions.jl/pull/1433) but it got stuck due to some ForwardDiff weirdness. I added some motivation in the PR there for why I became convinced that generally parameters and the return type of rand should be treated separately.