Remove eltype definitions, add deprecation warning, and update docu…

…mentation

docs/src/extends.md

@@ -4,11 +4,8 @@ Whereas this package already provides a large collection of common distributions
 
 Generally, you don't have to implement every API method listed in the documentation. This package provides a series of generic functions that turn a small number of internal methods into user-end API methods. What you need to do is to implement this small set of internal methods for your distributions.
 
-By default, `Discrete` sampleables have the support of type `Int` while `Continuous` sampleables have the support of type `Float64`. If this assumption does not hold for your new distribution or sampler, or its `ValueSupport` is neither `Discrete` nor `Continuous`, you should implement the `eltype` method in addition to the other methods listed below.
-
 **Note:** The methods that need to be implemented are different for distributions of different variate forms.
 
-
 ## Create a Sampler
 
 Unlike full-fledged distributions, a sampler, in general, only provides limited functionalities, mainly to support sampling.
@@ -18,60 +15,48 @@ Unlike full-fledged distributions, a sampler, in general, only provides limited
 To implement a univariate sampler, one can define a subtype (say `Spl`) of `Sampleable{Univariate,S}` (where `S` can be `Discrete` or `Continuous`), and provide a `rand` method, as
 
 ```julia
-function rand(rng::AbstractRNG, s::Spl)
+function Base.rand(rng::AbstractRNG, s::Spl)
     # ... generate a single sample from s
 end
 ```
 
-The package already implements a vectorized version of `rand!` and `rand` that repeatedly calls the scalar version to generate multiple samples; as wells as a one arg version that uses the default random number generator.
-
-### Multivariate Sampler
+The package already implements vectorized versions `rand!(rng::AbstractRNG, s::Spl, dims::Int...)` and `rand(rng::AbstractRNG, s::Spl, dims::Int...)` that repeatedly call the scalar version to generate multiple samples.
+Additionally, the package implements versions of these functions without the `rng::AbstractRNG` argument that use the default random number generator.
 
-To implement a multivariate sampler, one can define a subtype of `Sampleable{Multivariate,S}`, and provide both `length` and `_rand!` methods, as
+If there is a more efficient method to generate multiple samples, one should provide the following method
 
 ```julia
-Base.length(s::Spl) = ... # return the length of each sample
-
-function _rand!(rng::AbstractRNG, s::Spl, x::AbstractVector{T}) where T<:Real
-    # ... generate a single vector sample to x
+function Random.rand!(rng::AbstractRNG, s::Spl, x::AbstractArray{<:Real})
+    # ... generate multiple samples from s in x
 end
 ```
 
-This function can assume that the dimension of `x` is correct, and doesn't need to perform dimension checking.
+### Multivariate Sampler
 
-The package implements both `rand` and `rand!` as follows (which you don't need to implement in general):
+To implement a multivariate sampler, one can define a subtype of `Sampleable{Multivariate,S}`, and provide `length`, `rand`, and `rand!` methods, as
 
 ```julia
-function _rand!(rng::AbstractRNG, s::Sampleable{Multivariate}, A::DenseMatrix)
-    for i = 1:size(A,2)
-        _rand!(rng, s, view(A,:,i))
-    end
-    return A
-end
+Base.length(s::Spl) = ... # return the length of each sample
 
-function rand!(rng::AbstractRNG, s::Sampleable{Multivariate}, A::AbstractVector)
-    length(A) == length(s) ||
-        throw(DimensionMismatch("Output size inconsistent with sample length."))
-    _rand!(rng, s, A)
+function Base.rand(rng::AbstractRNG, s::Spl)
+    # ... generate a single vector sample from s
 end
 
-function rand!(rng::AbstractRNG, s::Sampleable{Multivariate}, A::DenseMatrix)
-    size(A,1) == length(s) ||
-        throw(DimensionMismatch("Output size inconsistent with sample length."))
-    _rand!(rng, s, A)
+@inline function Random.rand!(rng::AbstractRNG, s::Spl, x::AbstractVector{<:Real})
+    # `@inline` + `@boundscheck` allows users to skip bound checks by calling `@inbounds rand!(...)`
+    # Ref https://docs.julialang.org/en/v1/devdocs/boundscheck/#Eliding-bounds-checks
+    @boundscheck # ... check size (and possibly indices) of `x`
+    # ... generate a single vector sample from s in x
 end
-
-rand(rng::AbstractRNG, s::Sampleable{Multivariate,S}) where {S<:ValueSupport} =
-    _rand!(rng, s, Vector{eltype(S)}(length(s)))
-
-rand(rng::AbstractRNG, s::Sampleable{Multivariate,S}, n::Int) where {S<:ValueSupport} =
-    _rand!(rng, s, Matrix{eltype(S)}(length(s), n))
 ```
 
 If there is a more efficient method to generate multiple vector samples in a batch, one should provide the following method
 
 ```julia
-function _rand!(rng::AbstractRNG, s::Spl, A::DenseMatrix{T}) where T<:Real
+@inline function Random.rand!(rng::AbstractRNG, s::Spl, A::AbstractMatrix{<:Real})
+    # `@inline` + `@boundscheck` allows users to skip bound checks by calling `@inbounds rand!(...)`
+    # Ref https://docs.julialang.org/en/v1/devdocs/boundscheck/#Eliding-bounds-checks
+    @boundscheck # ... check size (and possibly indices) of `x`
     # ... generate multiple vector samples in batch
 end
 ```
@@ -80,17 +65,22 @@ Remember that each *column* of A is a sample.
 
 ### Matrix-variate Sampler
 
-To implement a multivariate sampler, one can define a subtype of `Sampleable{Multivariate,S}`, and provide both `size` and `_rand!` methods, as
+To implement a multivariate sampler, one can define a subtype of `Sampleable{Multivariate,S}`, and provide `size`, `rand`, and `rand!` methods, as
 
 ```julia
 Base.size(s::Spl) = ... # the size of each matrix sample
 
-function _rand!(rng::AbstractRNG, s::Spl, x::DenseMatrix{T}) where T<:Real
-    # ... generate a single matrix sample to x
+function Base.rand(rng::AbstractRNG, s::Spl)
+    # ... generate a single matrix sample from s
 end
-```
 
-Note that you can assume `x` has correct dimensions in `_rand!` and don't have to perform dimension checking, the generic `rand` and `rand!` will do dimension checking and array allocation for you.
+@inline function Random.rand!(rng::AbstractRNG, s::Spl, x::AbstractMatrix{<:Real})
+    # `@inline` + `@boundscheck` allows users to skip bound checks by calling `@inbounds rand!(...)`
+    # Ref https://docs.julialang.org/en/v1/devdocs/boundscheck/#Eliding-bounds-checks
+    @boundscheck # ... check size (and possibly indices) of `x`
+    # ... generate a single matrix sample from s in x
+end
+```
 
 ## Create a Distribution
 
@@ -106,7 +96,7 @@ A univariate distribution type should be defined as a subtype of `DiscreteUnivar
 
 The following methods need to be implemented for each univariate distribution type:
 
-- [`rand(::AbstractRNG, d::UnivariateDistribution)`](@ref)
+- [`Base.rand(::AbstractRNG, d::UnivariateDistribution)`](@ref)
 - [`sampler(d::Distribution)`](@ref)
 - [`logpdf(d::UnivariateDistribution, x::Real)`](@ref)
 - [`cdf(d::UnivariateDistribution, x::Real)`](@ref)
@@ -138,8 +128,8 @@ The following methods need to be implemented for each multivariate distribution
 
 - [`length(d::MultivariateDistribution)`](@ref)
 - [`sampler(d::Distribution)`](@ref)
-- [`eltype(d::Distribution)`](@ref)
-- [`Distributions._rand!(::AbstractRNG, d::MultivariateDistribution, x::AbstractArray)`](@ref)
+- [`Base.rand(::AbstractRNG, d::MultivariateDistribution)`](@ref)
+- [`Random.rand!(::AbstractRNG, d::MultivariateDistribution, x::AbstractVector{<:Real})`](@ref)
 - [`Distributions._logpdf(d::MultivariateDistribution, x::AbstractArray)`](@ref)
 
 Note that if there exist faster methods for batch evaluation, one should override `_logpdf!` and `_pdf!`.
@@ -161,6 +151,7 @@ A matrix-variate distribution type should be defined as a subtype of `DiscreteMa
 The following methods need to be implemented for each matrix-variate distribution type:
 
 - [`size(d::MatrixDistribution)`](@ref)
-- [`Distributions._rand!(rng::AbstractRNG, d::MatrixDistribution, A::AbstractMatrix)`](@ref)
+- [`Base.rand(rng::AbstractRNG, d::MatrixDistribution)`](@ref)
+- [`Random.rand!(rng::AbstractRNG, d::MatrixDistribution, A::AbstractMatrix{<:Real})`](@ref)
 - [`sampler(d::MatrixDistribution)`](@ref)
 - [`Distributions._logpdf(d::MatrixDistribution, x::AbstractArray)`](@ref)

docs/src/multivariate.md

@@ -18,7 +18,6 @@ The methods listed below are implemented for each multivariate distribution, whi
 ```@docs
 length(::MultivariateDistribution)
 size(::MultivariateDistribution)
-eltype(::Type{MultivariateDistribution})
 mean(::MultivariateDistribution)
 var(::MultivariateDistribution)
 cov(::MultivariateDistribution)

docs/src/types.md

@@ -57,7 +57,6 @@ The basic functionalities that a sampleable object provides are to *retrieve inf
 length(::Sampleable)
 size(::Sampleable)
 nsamples(::Type{Sampleable}, ::Any)
-eltype(::Type{Sampleable})
 rand(::AbstractRNG, ::Sampleable)
 rand!(::AbstractRNG, ::Sampleable, ::AbstractArray)
 ```

src/censored.jl

@@ -112,8 +112,6 @@ function partype(d::Censored{<:UnivariateDistribution,<:ValueSupport,T}) where {
     return promote_type(partype(d.uncensored), T)
 end
 
-Base.eltype(::Type{<:Censored{D,S,T}}) where {D,S,T} = promote_type(T, eltype(D))
-
 #### Range and Support
 
 isupperbounded(d::LeftCensored) = isupperbounded(d.uncensored)

src/cholesky/lkjcholesky.jl

@@ -82,8 +82,6 @@ end
 #  Properties
 #  -----------------------------------------------------------------------------
 
-Base.eltype(::Type{LKJCholesky{T}}) where {T} = T
-
 function Base.size(d::LKJCholesky)
     p = d.d
     return (p, p)

src/common.jl

@@ -94,14 +94,21 @@ Base.size(s::Sampleable{Univariate}) = ()
 Base.size(s::Sampleable{Multivariate}) = (length(s),)
 
 """
-    eltype(::Type{Sampleable})
+    eltype(::Type{S}) where {S<:Distributions.Sampleable}
+
+The default element type of a sample from a sampler of type `S`.
+
+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!`.
+
+!!! warn
+    This method is deprecated and will be removed in an upcoming breaking release.
 
-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
+function Base.eltype(::Type{S}) where {S<:Sampleable}
+    Base.depwarn("`eltype(::Type{<:Distributions.Sampleable})` is deprecated and will be removed", :eltype)
+    return eltype(Base.promote_op(rand, S))
+end
 
 """
     nsamples(s::Sampleable)

src/multivariate/dirichlet.jl

@@ -50,8 +50,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))

src/multivariate/jointorderstatistics.jl

@@ -88,8 +88,6 @@ maximum(d::JointOrderStatistics) = Fill(maximum(d.dist), length(d))
 
 params(d::JointOrderStatistics) = tuple(params(d.dist)..., d.n, d.ranks)
 partype(d::JointOrderStatistics) = partype(d.dist)
-Base.eltype(::Type{<:JointOrderStatistics{D}}) where {D} = Base.eltype(D)
-Base.eltype(d::JointOrderStatistics) = eltype(d.dist)
 
 function logpdf(d::JointOrderStatistics, x::AbstractVector{<:Real})
     n = d.n

src/multivariate/mvlogitnormal.jl

@@ -52,8 +52,6 @@ canonform(d::MvLogitNormal{<:MvNormal}) = MvLogitNormal(canonform(d.normal))
 # Properties
 
 length(d::MvLogitNormal) = length(d.normal) + 1
-Base.eltype(::Type{<:MvLogitNormal{D}}) where {D} = eltype(D)
-Base.eltype(d::MvLogitNormal) = eltype(d.normal)
 params(d::MvLogitNormal) = params(d.normal)
 @inline partype(d::MvLogitNormal) = partype(d.normal)
 

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))

src/multivariate/mvnormal.jl

@@ -223,8 +223,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.Σ))

src/multivariate/mvnormalcanon.jl

@@ -154,7 +154,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))

src/multivariate/mvtdist.jl

@@ -101,7 +101,6 @@ 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
 
 # For entropy calculations see "Multivariate t Distributions and their Applications", S. Kotz & S. Nadarajah
 function entropy(d::GenericMvTDist)

src/multivariate/product.jl

@@ -31,10 +31,6 @@ function Product(v::V) where {S<:ValueSupport,T<:UnivariateDistribution{S},V<:Ab
 end
 
 length(d::Product) = length(d.v)
-function Base.eltype(::Type{<:Product{S,T}}) where {S<:ValueSupport,
-                                                    T<:UnivariateDistribution{S}}
-    return eltype(T)
-end
 
 _rand!(rng::AbstractRNG, d::Product, x::AbstractVector{<:Real}) =
     map!(Base.Fix1(rand, rng), x, d.v)

src/product.jl

@@ -70,10 +70,6 @@ const ArrayOfUnivariateDistribution{N,D,S<:ValueSupport,T} = ProductDistribution
 const FillArrayOfUnivariateDistribution{N,D<:Fill{<:Any,N},S<:ValueSupport,T} = ProductDistribution{N,0,D,S,T}
 
 ## General definitions
-function Base.eltype(::Type{<:ProductDistribution{<:Any,<:Any,<:Any,<:ValueSupport,T}}) where {T}
-    return T
-end
-
 size(d::ProductDistribution) = d.size
 
 mean(d::ProductDistribution) = reshape(mapreduce(vec ∘ mean, vcat, d.dists), size(d))

src/reshaped.jl

@@ -24,7 +24,6 @@ function _reshape_check_dims(dist::Distribution{<:ArrayLikeVariate}, dims::Dims)
 end
 
 Base.size(d::ReshapedDistribution) = d.dims
-Base.eltype(::Type{ReshapedDistribution{<:Any,<:ValueSupport,D}}) where {D} = eltype(D)
 
 partype(d::ReshapedDistribution) = partype(d.dist)
 params(d::ReshapedDistribution) = (d.dist, d.dims)

src/truncate.jl

@@ -126,9 +126,6 @@ end
 params(d::Truncated) = tuple(params(d.untruncated)..., d.lower, d.upper)
 partype(d::Truncated{<:UnivariateDistribution,<:ValueSupport,T}) where {T<:Real} = promote_type(partype(d.untruncated), T)
 
-Base.eltype(::Type{<:Truncated{D}}) where {D<:UnivariateDistribution} = eltype(D)
-Base.eltype(d::Truncated) = eltype(d.untruncated)
-
 ### range and support
 
 islowerbounded(d::RightTruncated) = islowerbounded(d.untruncated)

src/univariate/continuous/gumbel.jl

@@ -41,8 +41,6 @@ Gumbel(μ::Real=0.0) = Gumbel(μ, one(μ); 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(θ))

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.μ

src/univariate/discrete/bernoulli.jl

@@ -40,8 +40,6 @@ Bernoulli() = Bernoulli{Float64}(0.5)
 
 @distr_support Bernoulli false true
 
-Base.eltype(::Type{<:Bernoulli}) = Bool
-
 #### Conversions
 convert(::Type{Bernoulli{T}}, p::Real) where {T<:Real} = Bernoulli(T(p))
 Base.convert(::Type{Bernoulli{T}}, d::Bernoulli) where {T<:Real} = Bernoulli{T}(T(d.p))

src/univariate/discrete/bernoullilogit.jl

@@ -24,8 +24,6 @@ BernoulliLogit() = BernoulliLogit(0.0)
 
 @distr_support BernoulliLogit false true
 
-Base.eltype(::Type{<:BernoulliLogit}) = Bool
-
 #### Conversions
 Base.convert(::Type{BernoulliLogit{T}}, d::BernoulliLogit) where {T<:Real} = BernoulliLogit{T}(T(d.logitp))
 Base.convert(::Type{BernoulliLogit{T}}, d::BernoulliLogit{T}) where {T<:Real} = d

src/univariate/discrete/dirac.jl

@@ -22,8 +22,6 @@ struct Dirac{T} <: DiscreteUnivariateDistribution
     value::T
 end
 
-Base.eltype(::Type{Dirac{T}}) where {T} = T
-
 insupport(d::Dirac, x::Real) = x == d.value
 minimum(d::Dirac) = d.value
 maximum(d::Dirac) = d.value

src/univariate/discrete/discretenonparametric.jl

@@ -39,8 +39,6 @@ DiscreteNonParametric(vs::AbstractVector{T}, ps::AbstractVector{P}; check_args::
         T<:Real,P<:Real} =
     DiscreteNonParametric{T,P,typeof(vs),typeof(ps)}(vs, ps; check_args=check_args)
 
-Base.eltype(::Type{<:DiscreteNonParametric{T}}) where T = T
-
 # Conversion
 convert(::Type{DiscreteNonParametric{T,P,Ts,Ps}}, d::DiscreteNonParametric) where {T,P,Ts,Ps} =
     DiscreteNonParametric{T,P,Ts,Ps}(convert(Ts, support(d)), convert(Ps, probs(d)), check_args=false)

src/univariate/locationscale.jl

@@ -71,8 +71,6 @@ end
 const ContinuousAffineDistribution{T<:Real,D<:ContinuousUnivariateDistribution} = AffineDistribution{T,Continuous,D}
 const DiscreteAffineDistribution{T<:Real,D<:DiscreteUnivariateDistribution} = AffineDistribution{T,Discrete,D}
 
-Base.eltype(::Type{<:AffineDistribution{T}}) where T = T
-
 minimum(d::AffineDistribution) =
     d.σ > 0 ? d.μ + d.σ * minimum(d.ρ) : d.μ + d.σ * maximum(d.ρ)
 maximum(d::AffineDistribution) =

src/univariate/orderstatistic.jl

@@ -58,8 +58,6 @@ insupport(d::OrderStatistic, x::Real) = insupport(d.dist, x)
 
 params(d::OrderStatistic) = tuple(params(d.dist)..., d.n, d.rank)
 partype(d::OrderStatistic) = partype(d.dist)
-Base.eltype(::Type{<:OrderStatistic{D}}) where {D} = Base.eltype(D)
-Base.eltype(d::OrderStatistic) = eltype(d.dist)
 
 # distribution of the ith order statistic from an IID uniform distribution, with CDF Uᵢₙ(x)
 function _uniform_orderstatistic(d::OrderStatistic)