From 6cb2b58856ce297225256eb441e69be6bf7399b0 Mon Sep 17 00:00:00 2001
From: David Widmann <dev+github@devmotion.de>
Date: Mon, 22 Nov 2021 15:29:55 +0100
Subject: [PATCH] Fix `eltype` and allow to specify type in `rand`

---
 src/cholesky/lkjcholesky.jl                   |  9 ++--
 src/common.jl                                 |  7 +--
 src/genericrand.jl                            | 47 ++++++++++++-------
 src/matrixvariates.jl                         | 14 +++---
 src/mixtures/mixturemodel.jl                  | 10 ++--
 src/mixtures/unigmm.jl                        | 15 +++---
 src/multivariate/dirichlet.jl                 |  2 -
 src/multivariate/mvlognormal.jl               |  2 -
 src/multivariate/mvnormal.jl                  |  6 +--
 src/multivariate/mvnormalcanon.jl             |  1 -
 src/multivariate/mvtdist.jl                   |  3 +-
 src/multivariates.jl                          | 22 ++++-----
 src/samplers/aliastable.jl                    |  4 +-
 src/samplers/binomial.jl                      | 33 ++++++-------
 src/samplers/categorical.jl                   |  6 +--
 src/samplers/discretenonparametric.jl         |  6 ++-
 src/samplers/exponential.jl                   |  6 ++-
 src/samplers/gamma.jl                         | 14 +++---
 src/samplers/poisson.jl                       | 24 +++++-----
 src/samplers/poissonbinomial.jl               |  4 +-
 src/samplers/vonmises.jl                      |  2 +-
 src/truncate.jl                               | 10 ++--
 src/truncated/normal.jl                       |  6 +--
 src/univariate/continuous/exponential.jl      |  2 +-
 src/univariate/continuous/gumbel.jl           |  6 +--
 src/univariate/continuous/lognormal.jl        |  4 +-
 src/univariate/continuous/noncentralf.jl      | 14 ++----
 src/univariate/continuous/normal.jl           |  6 +--
 src/univariate/continuous/tdist.jl            |  8 +++-
 src/univariate/continuous/uniform.jl          |  2 -
 src/univariate/discrete/dirac.jl              |  9 +++-
 .../discrete/discretenonparametric.jl         |  4 +-
 src/univariate/discrete/negativebinomial.jl   |  4 +-
 src/univariate/discrete/poisson.jl            |  2 +-
 src/univariate/locationscale.jl               | 18 +++++--
 src/univariates.jl                            | 20 ++++----
 36 files changed, 190 insertions(+), 162 deletions(-)

diff --git a/src/cholesky/lkjcholesky.jl b/src/cholesky/lkjcholesky.jl
index 12beecfd8..2f2a6433b 100644
--- a/src/cholesky/lkjcholesky.jl
+++ b/src/cholesky/lkjcholesky.jl
@@ -79,8 +79,6 @@ end
 #  Properties
 #  -----------------------------------------------------------------------------
 
-Base.eltype(::Type{LKJCholesky{T}}) where {T} = T
-
 function Base.size(d::LKJCholesky)
     p = d.d
     return (p, p)
@@ -150,15 +148,14 @@ end
 #  Sampling
 #  -----------------------------------------------------------------------------
 
-function Base.rand(rng::AbstractRNG, d::LKJCholesky)
-    factors = Matrix{eltype(d)}(undef, size(d))
+function Base.rand(rng::AbstractRNG, ::Type{T}, d::LKJCholesky) where {T}
+    factors = Matrix{T}(undef, size(d))
     R = LinearAlgebra.Cholesky(factors, d.uplo, 0)
     return _lkj_cholesky_onion_sampler!(rng, d, R)
 end
-function Base.rand(rng::AbstractRNG, d::LKJCholesky, dims::Dims)
+function Base.rand(rng::AbstractRNG, ::Type{T}, d::LKJCholesky, dims::Dims) where {T}
     p = d.d
     uplo = d.uplo
-    T = eltype(d)
     TM = Matrix{T}
     Rs = Array{LinearAlgebra.Cholesky{T,TM}}(undef, dims)
     for i in eachindex(Rs)
diff --git a/src/common.jl b/src/common.jl
index 243e1a591..9f5c4ed00 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -64,9 +64,10 @@ Base.size(s::Sampleable{Multivariate}) = (length(s),)
 """
     eltype(::Type{Sampleable})
 
-The default element type of a sample. This is the type of elements of the samples generated
-by the `rand` method. However, one can provide an array of different element types to
-store the samples using `rand!`.
+The default element type of a sample.
+
+This is the type of elements of the samples generated by the `rand` method. However, one can
+provide an array of different element types to store the samples using `rand!`.
 """
 Base.eltype(::Type{<:Sampleable{F,Discrete}}) where {F} = Int
 Base.eltype(::Type{<:Sampleable{F,Continuous}}) where {F} = Float64
diff --git a/src/genericrand.jl b/src/genericrand.jl
index f54c2748e..ce013cb32 100644
--- a/src/genericrand.jl
+++ b/src/genericrand.jl
@@ -1,30 +1,42 @@
 ### Generic rand methods
 
 """
-    rand([rng::AbstractRNG,] s::Sampleable)
+    rand(rng::AbstractRNG=GLOBAL_RNG, ::Type{T}=eltype(s), s::Sampleable)
 
-Generate one sample for `s`.
+Generate one sample for `s` of elment type `T`.
+"""
+rand(s::Sampleable) = rand(eltype(s), s)
+rand(::Type{T}, s::Sampleable) where {T} = rand(GLOBAL_RNG, T, s)
 
-    rand([rng::AbstractRNG,] s::Sampleable, n::Int)
+"""
+    rand(rng::AbstractRNG=GLOBAL_RNG, ::Type{T}=eltype(s), s::Sampleable, n::Int)
 
-Generate `n` samples from `s`. The form of the returned object depends on the variate form of `s`:
+Generate `n` samples from `s` of element type `T`.
 
+The form of the returned object depends on the variate form of `s`:
 - When `s` is univariate, it returns a vector of length `n`.
 - When `s` is multivariate, it returns a matrix with `n` columns.
 - When `s` is matrix-variate, it returns an array, where each element is a sample matrix.
+"""
+rand(s::Sampleable, n::Int) = rand(eltype(s), s, n)
+rand(::Type{T}, s::Sampleable, n::Int) where {T} = rand(GLOBAL_RNG, T, s, n)
+rand(rng::AbstractRNG, ::Type{T}, s::Sampleable, n::Int) where {T} = rand(rng, T, s, (n,))
 
-    rand([rng::AbstractRNG,] s::Sampleable, dim1::Int, dim2::Int...)
-    rand([rng::AbstractRNG,] s::Sampleable, dims::Dims)
+"""
+    rand(rng::AbstractRNG=GLOBAL_RNG, ::Type{T}=eltype(s), s::Sampleable, dims::Int...)
+    rand(rng::AbstractRNG=GLOBAL_RNG, ::Type{T}=eltype(s), s::Sampleable, dims::Dims)
 
-Generate an array of samples from `s` whose shape is determined by the given
+Generate an array of samples from `s` of element type `T` whose shape is determined by the given
 dimensions.
 """
-rand(s::Sampleable) = rand(GLOBAL_RNG, s)
-rand(s::Sampleable, dims::Dims) = rand(GLOBAL_RNG, s, dims)
-rand(s::Sampleable, dim1::Int, moredims::Int...) =
-    rand(GLOBAL_RNG, s, (dim1, moredims...))
-rand(rng::AbstractRNG, s::Sampleable, dim1::Int, moredims::Int...) =
-    rand(rng, s, (dim1, moredims...))
+rand(s::Sampleable, dims::Dims) = rand(eltype(s), s, dims)
+function rand(s::Sampleable, dims1::Int, dims2::Int, dims::Int...)
+    return rand(eltype(s), s, dims1, dims2, dims...)
+end
+function rand(::Type{T}, s::Sampleable, dims1::Int, dims2::Int, dims::Int...) where {T}
+    return rand(T, s, (dims1, dims2, dims...))
+end
+rand(::Type{T}, s::Sampleable, dims::Dims) where {T} = rand(GLOBAL_RNG, T, s, dims)
 
 """
     rand!([rng::AbstractRNG,] s::Sampleable, A::AbstractArray)
@@ -46,13 +58,12 @@ rand!(s::Sampleable, X::AbstractArray) = rand!(GLOBAL_RNG, s, X)
 rand!(rng::AbstractRNG, s::Sampleable, X::AbstractArray) = _rand!(rng, s, X)
 
 """
-    sampler(d::Distribution) -> Sampleable
-    sampler(s::Sampleable) -> s
+    sampler(s::Sampleable)
+
+Return a sampler that is used for batch sampling.
 
 Samplers can often rely on pre-computed quantities (that are not parameters
-themselves) to improve efficiency. If such a sampler exists, it can be provided
-with this `sampler` method, which would be used for batch sampling.
-The general fallback is `sampler(d::Distribution) = d`.
+themselves) to improve efficiency. The general fallback is `sampler(s) = s`.
 """
 sampler(s::Sampleable) = s
 
diff --git a/src/matrixvariates.jl b/src/matrixvariates.jl
index c7ac4ed4d..adfb69a96 100644
--- a/src/matrixvariates.jl
+++ b/src/matrixvariates.jl
@@ -128,16 +128,14 @@ function rand!(rng::AbstractRNG, s::Sampleable{Matrixvariate},
 end
 
 # multiple matrix-variates, must allocate array of arrays
-rand(rng::AbstractRNG, s::Sampleable{Matrixvariate}, dims::Dims) =
-    rand!(rng, s, Array{Matrix{eltype(s)}}(undef, dims), true)
-rand(rng::AbstractRNG, s::Sampleable{Matrixvariate,Continuous}, dims::Dims) =
-    rand!(rng, s, Array{Matrix{float(eltype(s))}}(undef, dims), true)
+function rand(rng::AbstractRNG, ::Type{T}, s::Sampleable{Matrixvariate}, dims::Dims) where {T}
+    return rand!(rng, s, Array{Matrix{T}}(undef, dims), true)
+end
 
 # single matrix-variate, must allocate one matrix
-rand(rng::AbstractRNG, s::Sampleable{Matrixvariate}) =
-    _rand!(rng, s, Matrix{eltype(s)}(undef, size(s)))
-rand(rng::AbstractRNG, s::Sampleable{Matrixvariate,Continuous}) =
-    _rand!(rng, s, Matrix{float(eltype(s))}(undef, size(s)))
+function rand(rng::AbstractRNG, ::Type{T}, s::Sampleable{Matrixvariate}) where {T}
+    return _rand!(rng, s, Matrix{T}(undef, size(s)))
+end
 
 # single matrix-variate with pre-allocated matrix
 function rand!(rng::AbstractRNG, s::Sampleable{Matrixvariate},
diff --git a/src/mixtures/mixturemodel.jl b/src/mixtures/mixturemodel.jl
index 7ced3ff1f..ca14a47a5 100644
--- a/src/mixtures/mixturemodel.jl
+++ b/src/mixtures/mixturemodel.jl
@@ -468,10 +468,12 @@ function MixtureSampler(d::MixtureModel{VF,VS}) where {VF,VS}
     MixtureSampler{VF,VS,eltype(csamplers)}(csamplers, psampler)
 end
 
-rand(rng::AbstractRNG, s::MixtureSampler{Univariate}) =
-    rand(rng, s.csamplers[rand(rng, s.psampler)])
-rand(rng::AbstractRNG, d::MixtureModel{Univariate}) =
-    rand(rng, component(d, rand(rng, d.prior)))
+function rand(rng::AbstractRNG, ::Type{T}, s::MixtureSampler{Univariate}) where {T}
+    return rand(rng, T, s.csamplers[rand(rng, s.psampler)])
+end
+function rand(rng::AbstractRNG, ::Type{T}, d::MixtureModel{Univariate}) where {T}
+    return rand(rng, T, component(d, rand(rng, d.prior)))
+end
 
 # multivariate mixture sampler for a vector
 _rand!(rng::AbstractRNG, s::MixtureSampler{Multivariate}, x::AbstractVector) =
diff --git a/src/mixtures/unigmm.jl b/src/mixtures/unigmm.jl
index f1d808ec1..77b20bb31 100644
--- a/src/mixtures/unigmm.jl
+++ b/src/mixtures/unigmm.jl
@@ -25,10 +25,10 @@ probs(d::UnivariateGMM) = probs(d.prior)
 
 mean(d::UnivariateGMM) = dot(d.means, probs(d))
 
-rand(d::UnivariateGMM) = (k = rand(d.prior); d.means[k] + randn() * d.stds[k])
-
-rand(rng::AbstractRNG, d::UnivariateGMM) =
-    (k = rand(rng, d.prior); d.means[k] + randn(rng) * d.stds[k])
+function rand(rng::AbstractRNG, ::Type{T}, d::UnivariateGMM) where {T}
+    k = rand(rng, d.prior)
+    return T(d.means[k]) + randn(rng, T) * T(d.stds[k])
+end
 
 params(d::UnivariateGMM) = (d.means, d.stds, d.prior)
 
@@ -38,6 +38,9 @@ struct UnivariateGMMSampler{VT1<:AbstractVector{<:Real},VT2<:AbstractVector{<:Re
     psampler::AliasTable
 end
 
-rand(rng::AbstractRNG, s::UnivariateGMMSampler) =
-    (k = rand(rng, s.psampler); s.means[k] + randn(rng) * s.stds[k])
+function rand(rng::AbstractRNG, ::Type{T}, s::UnivariateGMMSampler) where {T}
+    k = rand(rng, s.psampler)
+    return T(s.means[k]) + randn(rng, T) * T(s.stds[k])
+end
+
 sampler(d::UnivariateGMM) = UnivariateGMMSampler(d.means, d.stds, sampler(d.prior))
diff --git a/src/multivariate/dirichlet.jl b/src/multivariate/dirichlet.jl
index afac39e5d..befc06ccc 100644
--- a/src/multivariate/dirichlet.jl
+++ b/src/multivariate/dirichlet.jl
@@ -46,8 +46,6 @@ end
 
 length(d::DirichletCanon) = length(d.alpha)
 
-Base.eltype(::Type{<:Dirichlet{T}}) where {T} = T
-
 #### Conversions
 convert(::Type{Dirichlet{T}}, cf::DirichletCanon) where {T<:Real} =
     Dirichlet(convert(AbstractVector{T}, cf.alpha))
diff --git a/src/multivariate/mvlognormal.jl b/src/multivariate/mvlognormal.jl
index 25cddeae8..094977593 100644
--- a/src/multivariate/mvlognormal.jl
+++ b/src/multivariate/mvlognormal.jl
@@ -176,8 +176,6 @@ MvLogNormal(μ::AbstractVector,s::Real) = MvLogNormal(MvNormal(μ,s))
 MvLogNormal(σ::AbstractVector) = MvLogNormal(MvNormal(σ))
 MvLogNormal(d::Int,s::Real) = MvLogNormal(MvNormal(d,s))
 
-Base.eltype(::Type{<:MvLogNormal{T}}) where {T} = T
-
 ### Conversion
 function convert(::Type{MvLogNormal{T}}, d::MvLogNormal) where T<:Real
     MvLogNormal(convert(MvNormal{T}, d.normal))
diff --git a/src/multivariate/mvnormal.jl b/src/multivariate/mvnormal.jl
index 34f672b1a..9322a35e3 100644
--- a/src/multivariate/mvnormal.jl
+++ b/src/multivariate/mvnormal.jl
@@ -80,8 +80,8 @@ abstract type AbstractMvNormal <: ContinuousMultivariateDistribution end
 insupport(d::AbstractMvNormal, x::AbstractVector) =
     length(d) == length(x) && all(isfinite, x)
 
-minimum(d::AbstractMvNormal) = fill(eltype(d)(-Inf), length(d))
-maximum(d::AbstractMvNormal) = fill(eltype(d)(Inf), length(d))
+minimum(d::AbstractMvNormal) = fill(-Inf, length(d))
+maximum(d::AbstractMvNormal) = fill(Inf, length(d))
 mode(d::AbstractMvNormal) = mean(d)
 modes(d::AbstractMvNormal) = [mean(d)]
 
@@ -222,8 +222,6 @@ Base.@deprecate MvNormal(μ::AbstractVector{<:Real}, σ::Real) MvNormal(μ, σ^2
 Base.@deprecate MvNormal(σ::AbstractVector{<:Real}) MvNormal(LinearAlgebra.Diagonal(map(abs2, σ)))
 Base.@deprecate MvNormal(d::Int, σ::Real) MvNormal(LinearAlgebra.Diagonal(FillArrays.Fill(σ^2, d)))
 
-Base.eltype(::Type{<:MvNormal{T}}) where {T} = T
-
 ### Conversion
 function convert(::Type{MvNormal{T}}, d::MvNormal) where T<:Real
     MvNormal(convert(AbstractArray{T}, d.μ), convert(AbstractArray{T}, d.Σ))
diff --git a/src/multivariate/mvnormalcanon.jl b/src/multivariate/mvnormalcanon.jl
index 3dd2430b4..63b647c70 100644
--- a/src/multivariate/mvnormalcanon.jl
+++ b/src/multivariate/mvnormalcanon.jl
@@ -151,7 +151,6 @@ length(d::MvNormalCanon) = length(d.μ)
 mean(d::MvNormalCanon) = convert(Vector{eltype(d.μ)}, d.μ)
 params(d::MvNormalCanon) = (d.μ, d.h, d.J)
 @inline partype(d::MvNormalCanon{T}) where {T<:Real} = T
-Base.eltype(::Type{<:MvNormalCanon{T}}) where {T} = T
 
 var(d::MvNormalCanon) = diag(inv(d.J))
 cov(d::MvNormalCanon) = Matrix(inv(d.J))
diff --git a/src/multivariate/mvtdist.jl b/src/multivariate/mvtdist.jl
index edc9cec64..2f9f96694 100644
--- a/src/multivariate/mvtdist.jl
+++ b/src/multivariate/mvtdist.jl
@@ -98,8 +98,7 @@ invcov(d::GenericMvTDist) = d.df>2 ? ((d.df-2)/d.df)*Matrix(inv(d.Σ)) : NaN*one
 logdet_cov(d::GenericMvTDist) = d.df>2 ? logdet((d.df/(d.df-2))*d.Σ) : NaN
 
 params(d::GenericMvTDist) = (d.df, d.μ, d.Σ)
-@inline partype(d::GenericMvTDist{T}) where {T} = T
-Base.eltype(::Type{<:GenericMvTDist{T}}) where {T} = T
+@inline partype(::GenericMvTDist{T}) where {T} = T
 
 # For entropy calculations see "Multivariate t Distributions and their Applications", S. Kotz & S. Nadarajah
 function entropy(d::GenericMvTDist)
diff --git a/src/multivariates.jl b/src/multivariates.jl
index 0d2f214c9..a355f5dd0 100644
--- a/src/multivariates.jl
+++ b/src/multivariates.jl
@@ -66,21 +66,17 @@ function rand!(rng::AbstractRNG, s::Sampleable{Multivariate},
 end
 
 # multiple multivariate, must allocate matrix or array of vectors
-rand(s::Sampleable{Multivariate}, n::Int) = rand(GLOBAL_RNG, s, n)
-rand(rng::AbstractRNG, s::Sampleable{Multivariate}, n::Int) =
-    _rand!(rng, s, Matrix{eltype(s)}(undef, length(s), n))
-rand(rng::AbstractRNG, s::Sampleable{Multivariate,Continuous}, n::Int) =
-    _rand!(rng, s, Matrix{float(eltype(s))}(undef, length(s), n))
-rand(rng::AbstractRNG, s::Sampleable{Multivariate}, dims::Dims) =
-    rand!(rng, s, Array{Vector{eltype(s)}}(undef, dims), true)
-rand(rng::AbstractRNG, s::Sampleable{Multivariate,Continuous}, dims::Dims) =
-    rand!(rng, s, Array{Vector{float(eltype(s))}}(undef, dims), true)
+function rand(rng::AbstractRNG, ::Type{T}, s::Sampleable{Multivariate,Continuous}, n::Int) where {T}
+    return _rand!(rng, s, Matrix{T}(undef, length(s), n))
+end
+function rand(rng::AbstractRNG, ::Type{T}, s::Sampleable{Multivariate}, dims::Dims) where {T}
+    return rand!(rng, s, Array{Vector{T}}(undef, dims), true)
+end
 
 # single multivariate, must allocate vector
-rand(rng::AbstractRNG, s::Sampleable{Multivariate}) =
-    _rand!(rng, s, Vector{eltype(s)}(undef, length(s)))
-rand(rng::AbstractRNG, s::Sampleable{Multivariate,Continuous}) =
-    _rand!(rng, s, Vector{float(eltype(s))}(undef, length(s)))
+function rand(rng::AbstractRNG, ::Type{T}, s::Sampleable{Multivariate}) where {T}
+    return _rand!(rng, s, Vector{T}(undef, length(s)))
+end
 
 ## domain
 
diff --git a/src/samplers/aliastable.jl b/src/samplers/aliastable.jl
index ddebbce57..1c4e82a5f 100644
--- a/src/samplers/aliastable.jl
+++ b/src/samplers/aliastable.jl
@@ -13,11 +13,11 @@ function AliasTable(probs::AbstractVector)
     AliasTable(accp, alias)
 end
 
-function rand(rng::AbstractRNG, s::AliasTable)
+function rand(rng::AbstractRNG, ::Type{T}, s::AliasTable) where {T}
     i = rand(rng, 1:length(s.alias)) % Int
     u = rand(rng)
     @inbounds r = u < s.accept[i] ? i : s.alias[i]
-    r
+    return T(r)
 end
 
 show(io::IO, s::AliasTable) = @printf(io, "AliasTable with %d entries", ncategories(s))
diff --git a/src/samplers/binomial.jl b/src/samplers/binomial.jl
index 340a07424..6450b8cbc 100644
--- a/src/samplers/binomial.jl
+++ b/src/samplers/binomial.jl
@@ -44,23 +44,23 @@ function BinomialGeomSampler(n::Int, prob::Float64)
     BinomialGeomSampler(comp, n, scale)
 end
 
-function rand(rng::AbstractRNG, s::BinomialGeomSampler)
-    y = 0
-    x = 0
-    n = s.n
+function rand(rng::AbstractRNG, ::Type{T}, s::BinomialGeomSampler) where {T}
+    y = zero(T)
+    x = zero(T)
+    n = T(s.n)
     while true
         er = randexp(rng)
         v = er * s.scale
         if v > n  # in case when v is very large or infinity
             break
         end
-        y += ceil(Int,v)
+        y += T(ceil(v))
         if y > n
             break
         end
         x += 1
     end
-    (s.comp ? s.n - x : x)::Int
+    return s.comp ? n - x : x
 end
 
 
@@ -129,14 +129,14 @@ function BinomialTPESampler(n::Int, prob::Float64)
                        xM,xL,xR,c,λL,λR)
 end
 
-function rand(rng::AbstractRNG, s::BinomialTPESampler)
-    y = 0
+function rand(rng::AbstractRNG, ::Type{T}, s::BinomialTPESampler) where {T}
+    y = zero(T)
     while true
         # Step 1
         u = s.p4*rand(rng)
         v = rand(rng)
         if u <= s.p1
-            y = floor(Int,s.xM-s.p1*v+u)
+            y = T(floor(s.xM-s.p1*v+u))
             # Goto 6
             break
         elseif u <= s.p2 # Step 2
@@ -146,10 +146,10 @@ function rand(rng::AbstractRNG, s::BinomialTPESampler)
                 # Goto 1
                 continue
             end
-            y = floor(Int,x)
+            y = T(floor(x))
             # Goto 5
         elseif u <= s.p3 # Step 3
-            y = floor(Int,s.xL + log(v)/s.λL)
+            y = T(floor(s.xL + log(v)/s.λL))
             if y < 0
                 # Goto 1
                 continue
@@ -157,7 +157,7 @@ function rand(rng::AbstractRNG, s::BinomialTPESampler)
             v *= (u-s.p2)*s.λL
             # Goto 5
         else # Step 4
-            y = floor(Int,s.xR-log(v)/s.λR)
+            y = T(floor(s.xR-log(v)/s.λR))
             if y > s.n
                 # Goto 1
                 continue
@@ -219,7 +219,7 @@ function rand(rng::AbstractRNG, s::BinomialTPESampler)
         end
     end
     # 6
-    (s.comp ? s.n - y : y)::Int
+    return s.comp ? T(s.n) - y : y
 end
 
 
@@ -231,7 +231,7 @@ end
 
 BinomialAliasSampler(n::Int, p::Float64) = BinomialAliasSampler(AliasTable(binompvec(n, p)))
 
-rand(rng::AbstractRNG, s::BinomialAliasSampler) = rand(rng, s.table) - 1
+rand(rng::AbstractRNG, ::Type{T}, s::BinomialAliasSampler) where {T} = rand(rng, T, s.table) - 1
 
 
 # Integrated Polyalgorithm sampler that automatically chooses the proper one
@@ -260,5 +260,6 @@ end
 
 BinomialPolySampler(n::Real, p::Real) = BinomialPolySampler(round(Int, n), Float64(p))
 
-rand(rng::AbstractRNG, s::BinomialPolySampler) =
-    s.use_btpe ? rand(rng, s.btpe_sampler) : rand(rng, s.geom_sampler)
+function rand(rng::AbstractRNG, ::Type{T}, s::BinomialPolySampler) where {T}
+    return rand(rng, T, s.use_btpe ? s.btpe_sampler : s.geom_sampler)
+end
diff --git a/src/samplers/categorical.jl b/src/samplers/categorical.jl
index ae4f472d2..568616f34 100644
--- a/src/samplers/categorical.jl
+++ b/src/samplers/categorical.jl
@@ -15,14 +15,14 @@ CategoricalDirectSampler(p::Ts) where {T<:Real,Ts<:AbstractVector{T}} =
 
 ncategories(s::CategoricalDirectSampler) = length(s.prob)
 
-function rand(rng::AbstractRNG, s::CategoricalDirectSampler)
+function rand(rng::AbstractRNG, ::Type{T}, s::CategoricalDirectSampler) where {T}
     p = s.prob
     n = length(p)
     i = 1
     c = p[1]
-    u = rand(rng)
+    u = rand(rng, typeof(float(c)))
     while c < u && i < n
         c += p[i += 1]
     end
-    return i
+    return T(i)
 end
diff --git a/src/samplers/discretenonparametric.jl b/src/samplers/discretenonparametric.jl
index ce964a3ef..18ed2dfcf 100644
--- a/src/samplers/discretenonparametric.jl
+++ b/src/samplers/discretenonparametric.jl
@@ -19,5 +19,7 @@ DiscreteNonParametricSampler(support::S, probs::AbstractVector{<:Real}
     ) where {T<:Real,S<:AbstractVector{T}} =
     DiscreteNonParametricSampler{T,S}(support, probs)
 
-rand(rng::AbstractRNG, s::DiscreteNonParametricSampler) =
-    (@inbounds v = s.support[rand(rng, s.aliastable)]; v)
+function rand(rng::AbstractRNG, ::Type{T}, s::DiscreteNonParametricSampler) where {T}
+    @inbounds v = s.support[rand(rng, s.aliastable)]
+    return T(v)
+end
diff --git a/src/samplers/exponential.jl b/src/samplers/exponential.jl
index ccb335ebc..74d896761 100644
--- a/src/samplers/exponential.jl
+++ b/src/samplers/exponential.jl
@@ -2,10 +2,12 @@ struct ExponentialSampler <: Sampleable{Univariate,Continuous}
     scale::Float64
 end
 
-rand(rng::AbstractRNG, s::ExponentialSampler) = s.scale * randexp(rng)
+rand(rng::AbstractRNG, ::Type{T}, s::ExponentialSampler) where {T} = T(s.scale) * randexp(rng, T)
 
 struct ExponentialLogUSampler <: Sampleable{Univariate,Continuous}
     scale::Float64
 end
 
-rand(rng::AbstractRNG, s::ExponentialLogUSampler) = -s.scale * log(rand(rng))
+function rand(rng::AbstractRNG, ::Type{T}, s::ExponentialLogUSampler) where {T}
+    return -T(s.scale) * log(rand(rng, T))
+end
diff --git a/src/samplers/gamma.jl b/src/samplers/gamma.jl
index 70d9d8b1a..2b119f794 100644
--- a/src/samplers/gamma.jl
+++ b/src/samplers/gamma.jl
@@ -75,7 +75,7 @@ function calc_q(s::GammaGDSampler, t)
     end
 end
 
-function rand(rng::AbstractRNG, s::GammaGDSampler)
+function rand(rng::AbstractRNG, ::Type{Float64}, s::GammaGDSampler)
     # Step 2
     t = randn(rng)
     x = s.s + 0.5t
@@ -138,7 +138,7 @@ function GammaGSSampler(d::Gamma)
     GammaGSSampler(a, ia, b, scale(d))
 end
 
-function rand(rng::AbstractRNG, s::GammaGSSampler)
+function rand(rng::AbstractRNG, ::Type{Float64}, s::GammaGSSampler)
     while true
         # step 1
         p = s.b*rand(rng)
@@ -175,7 +175,7 @@ function GammaMTSampler(g::Gamma)
     GammaMTSampler(d, c, κ)
 end
 
-function rand(rng::AbstractRNG, s::GammaMTSampler)
+function rand(rng::AbstractRNG, ::Type{Float64}, s::GammaMTSampler)
     while true
         x = randn(rng)
         v = 1.0 + s.c * x
@@ -204,10 +204,10 @@ function GammaIPSampler(d::Gamma,::Type{S}) where S<:Sampleable
 end
 GammaIPSampler(d::Gamma) = GammaIPSampler(d,GammaMTSampler)
 
-function rand(rng::AbstractRNG, s::GammaIPSampler)
-    x = rand(rng, s.s)
-    e = randexp(rng)
-    x*exp(s.nia*e)
+function rand(rng::AbstractRNG, ::Type{T}, s::GammaIPSampler) where {T}
+    x = rand(rng, T, s.s)
+    e = randexp(rng, T)
+    return x * exp(T(s.nia) * e)
 end
 
 # function sampler(d::Gamma)
diff --git a/src/samplers/poisson.jl b/src/samplers/poisson.jl
index 64898ec2f..ef1ce6a2f 100644
--- a/src/samplers/poisson.jl
+++ b/src/samplers/poisson.jl
@@ -19,14 +19,14 @@ end
 
 PoissonCountSampler(d::Poisson) = PoissonCountSampler(rate(d))
 
-function rand(rng::AbstractRNG, s::PoissonCountSampler)
+function rand(rng::AbstractRNG, ::Type{T}, s::PoissonCountSampler) where {T}
     μ = s.μ
-    T = typeof(μ)
-    n = 0
-    c = randexp(rng, T)
+    n = zero(T)
+    C = typeof(float(μ))
+    c = randexp(rng, C)
     while c < μ
         n += 1
-        c += randexp(rng, T)
+        c += randexp(rng, C)
     end
     return n
 end
@@ -43,7 +43,7 @@ struct PoissonADSampler{T<:Real} <: Sampleable{Univariate,Discrete}
     μ::T
     s::T
     d::T
-    L::Int
+    L::T
 end
 
 PoissonADSampler(d::Poisson) = PoissonADSampler(rate(d))
@@ -51,23 +51,23 @@ PoissonADSampler(d::Poisson) = PoissonADSampler(rate(d))
 function PoissonADSampler(μ::Real)
     s = sqrt(μ)
     d = 6 * μ^2
-    L = floor(Int, μ - 1.1484)
+    L = floor(μ - 1.1484)
 
-    PoissonADSampler(promote(μ, s, d)..., L)
+    PoissonADSampler(promote(μ, s, d, L)...)
 end
 
-function rand(rng::AbstractRNG, sampler::PoissonADSampler)
+function rand(rng::AbstractRNG, ::Type{Ttype}, sampler::PoissonADSampler) where {Ttype}
     μ = sampler.μ
     s = sampler.s
     d = sampler.d
-    L = sampler.L
+    L = Ttype(sampler.L)
     μType = typeof(μ)
 
     # Step N
     G = μ + s * randn(rng, μType)
 
     if G >= zero(G)
-        K = floor(Int, G)
+        K = Ttype(floor(G))
         # Step I
         if K >= L
             return K
@@ -97,7 +97,7 @@ function rand(rng::AbstractRNG, sampler::PoissonADSampler)
             continue
         end
 
-        K = floor(Int, μ + s * T)
+        K = Ttype(floor(μ + s * T))
         px, py, fx, fy = procf(μ, K, s)
         c = 0.1069 / μ
 
diff --git a/src/samplers/poissonbinomial.jl b/src/samplers/poissonbinomial.jl
index 1c6b81502..e542d2fee 100644
--- a/src/samplers/poissonbinomial.jl
+++ b/src/samplers/poissonbinomial.jl
@@ -4,4 +4,6 @@ end
 
 PoissBinAliasSampler(d::PoissonBinomial) = PoissBinAliasSampler(AliasTable(d.pmf))
 
-rand(rng::AbstractRNG, s::PoissBinAliasSampler) = rand(rng, s.table) - 1
+function rand(rng::AbstractRNG, ::Type{T}, s::PoissBinAliasSampler) where {T}
+    return rand(rng, T, s.table) - 1
+end
diff --git a/src/samplers/vonmises.jl b/src/samplers/vonmises.jl
index b3d751b0f..b24d40415 100644
--- a/src/samplers/vonmises.jl
+++ b/src/samplers/vonmises.jl
@@ -15,7 +15,7 @@ end
 #     DJ Best & NI Fisher (1979). Efficient Simulation of the von Mises
 #     Distribution. Journal of the Royal Statistical Society. Series C
 #     (Applied Statistics), 28(2), 152-157.
-function rand(rng::AbstractRNG, s::VonMisesSampler)
+function rand(rng::AbstractRNG, ::Type{Float64}, s::VonMisesSampler)
     f = 0.0
     local x::Float64
     if s.κ > 700.0
diff --git a/src/truncate.jl b/src/truncate.jl
index ba76b7757..298e0019e 100644
--- a/src/truncate.jl
+++ b/src/truncate.jl
@@ -63,7 +63,7 @@ end
 
 params(d::Truncated) = tuple(params(d.untruncated)..., d.lower, d.upper)
 partype(d::Truncated) = partype(d.untruncated)
-Base.eltype(::Type{Truncated{D, S, T} } ) where {D, S, T} = T
+Base.eltype(::Type{<:Truncated{D}}) where {D} = eltype(D)
 
 ### range and support
 
@@ -137,17 +137,19 @@ end
 
 ## random number generation
 
-function _rand!(rng::AbstractRNG, d::Truncated)
+function rand(rng::AbstractRNG, ::Type{T}, d::Truncated) where {T}
     d0 = d.untruncated
     if d.tp > 0.25
         while true
-            r = _rand!(rng, d0)
+            r = rand(rng, T, d0)
             if d.lower <= r <= d.upper
                 return r
             end
         end
     else
-        return quantile(d0, d.lcdf + rand(rng) * d.tp)
+        lcdf = d.lcdf
+        tp = d.tp
+        return convert(T, quantile(d0, lcdf + rand(rng, typeof(lcdf)) * tp))
     end
 end
 
diff --git a/src/truncated/normal.jl b/src/truncated/normal.jl
index 5dc57a60d..7085b6bcc 100644
--- a/src/truncated/normal.jl
+++ b/src/truncated/normal.jl
@@ -137,7 +137,7 @@ end
 ## Use specialized sampler, as quantile-based method is inaccurate in
 ## tail regions of the Normal, issue #343
 
-function rand(rng::AbstractRNG, d::Truncated{Normal{T},Continuous}) where T <: Real
+function rand(rng::AbstractRNG, ::Type{T}, d::Truncated{<:Normal}) where {T}
     d0 = d.untruncated
     μ = mean(d0)
     σ = std(d0)
@@ -145,9 +145,9 @@ function rand(rng::AbstractRNG, d::Truncated{Normal{T},Continuous}) where T <: R
         a = (d.lower - μ) / σ
         b = (d.upper - μ) / σ
         z = randnt(rng, a, b, d.tp)
-        return μ + σ * z
+        return T(μ + σ * z)
     else
-        return clamp(μ, d.lower, d.upper)
+        return T(clamp(μ, d.lower, d.upper))
     end
 end
 
diff --git a/src/univariate/continuous/exponential.jl b/src/univariate/continuous/exponential.jl
index 5c395c962..f538a988a 100644
--- a/src/univariate/continuous/exponential.jl
+++ b/src/univariate/continuous/exponential.jl
@@ -98,7 +98,7 @@ cf(d::Exponential, t::Real) = 1/(1 - t * im * scale(d))
 
 
 #### Sampling
-rand(rng::AbstractRNG, d::Exponential) = xval(d, randexp(rng))
+rand(rng::AbstractRNG, ::Type{T}, d::Exponential) where {T} = T(d.θ) * randexp(rng, T)
 
 
 #### Fit model
diff --git a/src/univariate/continuous/gumbel.jl b/src/univariate/continuous/gumbel.jl
index 6be60d088..1392625b7 100644
--- a/src/univariate/continuous/gumbel.jl
+++ b/src/univariate/continuous/gumbel.jl
@@ -42,8 +42,6 @@ Gumbel() = Gumbel(0.0, 1.0, check_args=false)
 
 const DoubleExponential = Gumbel
 
-Base.eltype(::Type{Gumbel{T}}) where {T} = T
-
 #### Conversions
 
 convert(::Type{Gumbel{T}}, μ::S, θ::S) where {T <: Real, S <: Real} = Gumbel(T(μ), T(θ))
@@ -56,8 +54,8 @@ scale(d::Gumbel) = d.θ
 params(d::Gumbel) = (d.μ, d.θ)
 partype(::Gumbel{T}) where {T} = T
 
-function Base.rand(rng::Random.AbstractRNG, d::Gumbel)
-    return d.μ - d.θ * log(randexp(rng, float(eltype(d))))
+function Base.rand(rng::Random.AbstractRNG, ::Type{T}, d::Gumbel) where {T}
+    return T(d.μ) - T(d.θ) * log(randexp(rng, T))
 end
 function Random.rand!(rng::Random.AbstractRNG, d::Gumbel, x::AbstractArray)
     randexp!(rng, x)
diff --git a/src/univariate/continuous/lognormal.jl b/src/univariate/continuous/lognormal.jl
index 718851ea2..d9ea0cfe1 100644
--- a/src/univariate/continuous/lognormal.jl
+++ b/src/univariate/continuous/lognormal.jl
@@ -151,7 +151,9 @@ end
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::LogNormal) = exp(randn(rng) * d.σ + d.μ)
+function rand(rng::AbstractRNG, ::Type{T}, d::LogNormal) where {T}
+    return exp(randn(rng, T) * T(d.σ) + T(d.μ))
+end
 
 ## Fitting
 
diff --git a/src/univariate/continuous/noncentralf.jl b/src/univariate/continuous/noncentralf.jl
index 26fce4c76..b8fb3f39b 100644
--- a/src/univariate/continuous/noncentralf.jl
+++ b/src/univariate/continuous/noncentralf.jl
@@ -47,16 +47,8 @@ var(d::NoncentralF{T}) where {T<:Real} = d.ν2 > 4 ? 2d.ν2^2 *
 
 @_delegate_statsfuns NoncentralF nfdist ν1 ν2 λ
 
-function rand(rng::AbstractRNG, d::NoncentralF)
-    r1 = rand(rng, NoncentralChisq(d.ν1,d.λ)) / d.ν1
-    r2 = rand(rng, Chisq(d.ν2)) / d.ν2
-    r1 / r2
-end
-
-# TODO: remove RFunctions dependency once NoncentralChisq has its removed
-@rand_rdist(NoncentralF)
-function rand(d::NoncentralF)
-    r1 = rand(NoncentralChisq(d.ν1,d.λ)) / d.ν1
-    r2 = rand(Chisq(d.ν2)) / d.ν2
+function rand(rng::AbstractRNG, ::Type{T}, d::NoncentralF) where {T}
+    r1 = rand(rng, T, NoncentralChisq(d.ν1,d.λ)) / T(d.ν1)
+    r2 = rand(rng, T, Chisq(d.ν2)) / T(d.ν2)
     r1 / r2
 end
diff --git a/src/univariate/continuous/normal.jl b/src/univariate/continuous/normal.jl
index fe5feab84..c74b9a19b 100644
--- a/src/univariate/continuous/normal.jl
+++ b/src/univariate/continuous/normal.jl
@@ -60,8 +60,6 @@ params(d::Normal) = (d.μ, d.σ)
 location(d::Normal) = d.μ
 scale(d::Normal) = d.σ
 
-Base.eltype(::Type{Normal{T}}) where {T} = T
-
 #### Statistics
 
 mean(d::Normal) = d.μ
@@ -249,7 +247,9 @@ cf(d::Normal, t::Real) = exp(im * t * d.μ - d.σ^2 / 2 * t^2)
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::Normal{T}) where {T} = d.μ + d.σ * randn(rng, float(T))
+function rand(rng::AbstractRNG, ::Type{T}, d::Normal) where {T}
+    return T(d.μ) + T(d.σ) * randn(rng, T)
+end
 
 #### Fitting
 
diff --git a/src/univariate/continuous/tdist.jl b/src/univariate/continuous/tdist.jl
index e1ca3d700..889a910c6 100644
--- a/src/univariate/continuous/tdist.jl
+++ b/src/univariate/continuous/tdist.jl
@@ -78,7 +78,13 @@ end
 
 @_delegate_statsfuns TDist tdist ν
 
-rand(rng::AbstractRNG, d::TDist) = randn(rng) / ( isinf(d.ν) ? 1 : sqrt(rand(rng, Chisq(d.ν))/d.ν) )
+function rand(rng::AbstractRNG, ::Type{T}, d::TDist) where {T}
+    return if isinf(d.ν)
+        randn(rng, T)
+    else
+        randn(rng, T) / sqrt(rand(rng, T, Chisq(d.ν)) / T(d.ν))
+    end
+end
 
 function cf(d::TDist{T}, t::Real) where T <: Real
     isinf(d.ν) && return cf(Normal(zero(T), one(T)), t)
diff --git a/src/univariate/continuous/uniform.jl b/src/univariate/continuous/uniform.jl
index fc04a9546..61be4c9de 100644
--- a/src/univariate/continuous/uniform.jl
+++ b/src/univariate/continuous/uniform.jl
@@ -112,8 +112,6 @@ end
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::Uniform) = d.a + (d.b - d.a) * rand(rng)
-
 rand!(rng::AbstractRNG, d::Uniform, A::AbstractArray) =
     A .= quantile.(d, rand!(rng, A))
 
diff --git a/src/univariate/discrete/dirac.jl b/src/univariate/discrete/dirac.jl
index 7be41459f..83bf4c0bb 100644
--- a/src/univariate/discrete/dirac.jl
+++ b/src/univariate/discrete/dirac.jl
@@ -54,4 +54,11 @@ cf(d::Dirac, t) = cis(t * d.value)
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::Dirac) = d.value
+rand(::AbstractRNG, ::Type{T}, d::Dirac) where {T} = convert(T, d.value)
+function rand(::AbstractRNG, ::Type{T}, d::Dirac, dims::Dims) where {T}
+    return fill(convert(T, d.value), dims)
+end
+function rand!(::AbstractRNG, d::Dirac, x::AbstractArray)
+    fill!(x, d.value)
+    return x
+end
diff --git a/src/univariate/discrete/discretenonparametric.jl b/src/univariate/discrete/discretenonparametric.jl
index e54a0789f..dbfc87579 100644
--- a/src/univariate/discrete/discretenonparametric.jl
+++ b/src/univariate/discrete/discretenonparametric.jl
@@ -69,7 +69,7 @@ Base.isapprox(c1::D, c2::D) where D<:DiscreteNonParametric =
 
 # Sampling
 
-function rand(rng::AbstractRNG, d::DiscreteNonParametric)
+function rand(rng::AbstractRNG, ::Type{T}, d::DiscreteNonParametric) where {T}
     x = support(d)
     p = probs(d)
     n = length(p)
@@ -79,7 +79,7 @@ function rand(rng::AbstractRNG, d::DiscreteNonParametric)
     while cp <= draw && i < n
         @inbounds cp += p[i +=1]
     end
-    return x[i]
+    return convert(T, x[i])
 end
 
 sampler(d::DiscreteNonParametric) =
diff --git a/src/univariate/discrete/negativebinomial.jl b/src/univariate/discrete/negativebinomial.jl
index 2b0be1f58..0e59213cf 100644
--- a/src/univariate/discrete/negativebinomial.jl
+++ b/src/univariate/discrete/negativebinomial.jl
@@ -110,7 +110,9 @@ invlogcdf(d::NegativeBinomial, lq::Real) = convert(Int, nbinominvlogcdf(d.r, d.p
 invlogccdf(d::NegativeBinomial, lq::Real) = convert(Int, nbinominvlogccdf(d.r, d.p, lq))
 
 ## sampling
-rand(rng::AbstractRNG, d::NegativeBinomial) = rand(rng, Poisson(rand(rng, Gamma(d.r, (1 - d.p)/d.p))))
+function rand(rng::AbstractRNG, ::Type{T}, d::NegativeBinomial) where {T}
+    return rand(rng, T, Poisson(rand(rng, float(partype(d)), Gamma(d.r, (1 - d.p)/d.p))))
+end
 
 function mgf(d::NegativeBinomial, t::Real)
     r, p = params(d)
diff --git a/src/univariate/discrete/poisson.jl b/src/univariate/discrete/poisson.jl
index 75fa25b36..ccb8e7db1 100644
--- a/src/univariate/discrete/poisson.jl
+++ b/src/univariate/discrete/poisson.jl
@@ -143,4 +143,4 @@ function sampler(d::Poisson)
     end
 end
 
-rand(rng::AbstractRNG, d::Poisson) = rand(rng, sampler(d))
+rand(rng::AbstractRNG, ::Type{T}, d::Poisson) where {T} = rand(rng, T, sampler(d))
diff --git a/src/univariate/locationscale.jl b/src/univariate/locationscale.jl
index ba61c8d90..f8fb8416e 100644
--- a/src/univariate/locationscale.jl
+++ b/src/univariate/locationscale.jl
@@ -44,10 +44,9 @@ struct LocationScale{T<:Real, S<:ValueSupport, D<:UnivariateDistribution{S}} <:
 end
 
 function LocationScale(μ::T, σ::T, ρ::UnivariateDistribution; check_args=true) where {T<:Real}
-    _T = promote_type(eltype(ρ), T)
     D = typeof(ρ)
     S = value_support(D)
-    return LocationScale{_T,S,D}(_T(μ), _T(σ), ρ; check_args=check_args)
+    return LocationScale{T,S,D}(μ, σ, ρ; check_args=check_args)
 end
 
 LocationScale(μ::Real, σ::Real, ρ::UnivariateDistribution) = LocationScale(promote(μ, σ)..., ρ)
@@ -56,7 +55,18 @@ LocationScale(μ::Real, σ::Real, ρ::UnivariateDistribution) = LocationScale(pr
 const ContinuousLocationScale{T<:Real,D<:ContinuousUnivariateDistribution} = LocationScale{T,Continuous,D}
 const DiscreteLocationScale{T<:Real,D<:DiscreteUnivariateDistribution} = LocationScale{T,Discrete,D}
 
-Base.eltype(::Type{<:LocationScale{T}}) where T = T
+function Base.eltype(::Type{<:LocationScale{T,S,D}}) where {T,S,D<:ContinuousUnivariateDistribution}
+    return eltype(D)
+end
+function Base.eltype(::Type{<:LocationScale{T,S,D}}) where {T,S,D<:DiscreteUnivariateDistribution}
+    U = eltype(D)
+    return if float(U) === U
+        # no need to promote, already floating point number
+        U
+    else
+        promote(T, U)
+    end
+end
 
 minimum(d::LocationScale) = d.μ + d.σ * minimum(d.ρ)
 maximum(d::LocationScale) = d.μ + d.σ * maximum(d.ρ)
@@ -115,7 +125,7 @@ end
 
 quantile(d::LocationScale,q::Real) = d.μ + d.σ * quantile(d.ρ,q)
 
-rand(rng::AbstractRNG, d::LocationScale) = d.μ + d.σ * rand(rng, d.ρ)
+rand(rng::AbstractRNG, ::Type{T}, d::LocationScale) where {T} = convert(T, d.μ + d.σ * rand(rng, T, d.ρ))
 cf(d::LocationScale, t::Real) = cf(d.ρ,t*d.σ) * exp(1im*t*d.μ)
 gradlogpdf(d::ContinuousLocationScale, x::Real) = gradlogpdf(d.ρ,(x-d.μ)/d.σ) / d.σ
 
diff --git a/src/univariates.jl b/src/univariates.jl
index 52073863b..eba721c52 100644
--- a/src/univariates.jl
+++ b/src/univariates.jl
@@ -134,10 +134,9 @@ end
 ## sampling
 
 # multiple univariate, must allocate array
-rand(rng::AbstractRNG, s::Sampleable{Univariate}, dims::Dims) =
-    rand!(rng, s, Array{eltype(s)}(undef, dims))
-rand(rng::AbstractRNG, s::Sampleable{Univariate,Continuous}, dims::Dims) =
-    rand!(rng, s, Array{float(eltype(s))}(undef, dims))
+function rand(rng::AbstractRNG, ::Type{T}, s::Sampleable{Univariate}, dims::Dims) where {T}
+    return rand!(rng, s, Array{T}(undef, dims))
+end
 
 # multiple univariate with pre-allocated array
 # we use a function barrier since for some distributions `sampler(s)` is not type-stable:
@@ -147,18 +146,23 @@ function rand!(rng::AbstractRNG, s::Sampleable{Univariate}, A::AbstractArray)
 end
 
 function _rand_loops!(rng::AbstractRNG, sampler::Sampleable{Univariate}, A::AbstractArray)
+    T = eltype(A)
     for i in eachindex(A)
-        @inbounds A[i] = rand(rng, sampler)
+        @inbounds A[i] = rand(rng, T, sampler)
     end
     return A
 end
 
 """
-    rand(rng::AbstractRNG, d::UnivariateDistribution)
+    rand(rng::AbstractRNG=GLOBAL_RNG, ::Type{T}=eltype(d), d::UnivariateDistribution)
+
+Generate a scalar sample from `d` of element type `T`.
 
-Generate a scalar sample from `d`. The general fallback is `quantile(d, rand())`.
+The general fallback is `convert(T, quantile(d, rand(rng, float(T))))`.
 """
-rand(rng::AbstractRNG, d::UnivariateDistribution) = quantile(d, rand(rng))
+function rand(rng::AbstractRNG, ::Type{T}, d::UnivariateDistribution) where {T}
+    return convert(T, quantile(d, rand(rng, float(T))))
+end
 
 """
     rand!(rng::AbstractRNG, ::UnivariateDistribution, ::AbstractArray)