Generalize Product (#1391)

* Move `src/multivariate/product.jl`

* Generalize `Product` to `ProductDistribution`

* Add implementations for more general product distributions

* Unify and generalize `rand!`, `logpdf` and `pdf`

* Revert unrelated changes and fix tests

* Propagate `@inbounds`

* Remove unneeded implementation

* Fix typos

* Fix some dispatches

* More fixes

* Support tuple of distributions and mix of discrete + continuous

* Fix additional test errors

* Fix method ambiguity

* Fix `VonMisesFisherSampler`

* Fix mixture sampler

* Simplify multinomial sampler

* Fix `loglikelihood` for univariate distributions

* Add ReshapedDistribution

* Fix typo

* Revert some changes

* Update product.jl

* Remove duplicate `eachvariate`/`EachVariate`

* Reintroduce `Product`

* Improve type inference

* Add explanations of `ValueSupport`

* Fix typo

* Remove another breaking change

Author: devmotion

docs/src/types.md

@@ -33,16 +33,21 @@ The `VariateForm` sub-types defined in `Distributions.jl` are:
 
 ### ValueSupport
 
-```@doc
+```@docs
 Distributions.ValueSupport
 ```
 
 The `ValueSupport` sub-types defined in `Distributions.jl` are:
 
-**Type** | **Element type** | **Descriptions**
---- | --- | ---
-`Discrete` | `Int` | Samples take discrete values
-`Continuous` | `Float64` | Samples take continuous real values
+```@docs
+Distributions.Discrete
+Distributions.Continuous
+```
+
+**Type** | **Default element type** | **Description** | **Examples**
+--- | --- | --- | ---
+`Discrete` | `Int` | Samples take countably many values | $\{0,1,2,3\}$, $\mathbb{N}$
+`Continuous` | `Float64` | Samples take uncountably many values | $[0, 1]$, $\mathbb{R}$
 
 Multiple samples are often organized into an array, depending on the variate form.
 

src/Distributions.jl

@@ -145,7 +145,7 @@ export
     Pareto,
     PGeneralizedGaussian,
     SkewedExponentialPower,
-    Product,
+    Product, # deprecated
     Poisson,
     PoissonBinomial,
     QQPair,
@@ -293,6 +293,7 @@ include("cholesky/lkjcholesky.jl")
 include("samplers.jl")
 
 # others
+include("product.jl")
 include("reshaped.jl")
 include("truncate.jl")
 include("censored.jl")

src/common.jl

@@ -23,13 +23,42 @@ const Matrixvariate = ArrayLikeVariate{2}
 abstract type CholeskyVariate <: VariateForm end
 
 """
-`S <: ValueSupport` specifies the support of sample elements,
-either discrete or continuous.
+    ValueSupport
+
+Abstract type that specifies the support of elements of samples.
+
+It is either [`Discrete`](@ref) or [`Continuous`](@ref).
 """
 abstract type ValueSupport end
+
+"""
+    Discrete <: ValueSupport
+
+This type represents the support of a discrete random variable.
+
+It is countable. For instance, it can be a finite set or a countably infinite set such as
+the natural numbers.
+
+See also: [`Continuous`](@ref), [`ValueSupport`](@ref)
+"""
 struct Discrete   <: ValueSupport end
+
+"""
+    Continuous <: ValueSupport
+
+This types represents the support of a continuous random variable.
+
+It is uncountably infinite. For instance, it can be an interval on the real line.
+
+See also: [`Discrete`](@ref), [`ValueSupport`](@ref)
+"""
 struct Continuous <: ValueSupport end
 
+# promotions (e.g., in product distribution):
+# combination of discrete support (countable) and continuous support (uncountable) yields
+# continuous support (uncountable)
+Base.promote_rule(::Type{Continuous}, ::Type{Discrete}) = Continuous
+
 ## Sampleable
 
 """
@@ -42,7 +71,6 @@ Any `Sampleable` implements the `Base.rand` method.
 """
 abstract type Sampleable{F<:VariateForm,S<:ValueSupport} end
 
-
 variate_form(::Type{<:Sampleable{VF}}) where {VF} = VF
 value_support(::Type{<:Sampleable{<:VariateForm,VS}}) where {VS} = VS
 
@@ -142,10 +170,6 @@ const ContinuousMultivariateDistribution = Distribution{Multivariate,  Continuou
 const DiscreteMatrixDistribution         = Distribution{Matrixvariate, Discrete}
 const ContinuousMatrixDistribution       = Distribution{Matrixvariate, Continuous}
 
-variate_form(::Type{<:Distribution{VF}}) where {VF} = VF
-
-value_support(::Type{<:Distribution{VF,VS}}) where {VF,VS} = VS
-
 # allow broadcasting over distribution objects
 # to be decided: how to handle multivariate/matrixvariate distributions?
 Broadcast.broadcastable(d::UnivariateDistribution) = Ref(d)

src/multivariate/product.jl

@@ -1,4 +1,5 @@
-import Statistics: mean, var, cov
+# Deprecated product distribution
+# TODO: Remove in next breaking release
 
 """
     Product <: MultivariateDistribution
@@ -20,6 +21,10 @@ struct Product{
         V<:AbstractVector{T} where
         T<:UnivariateDistribution{S} where
         S<:ValueSupport
+        Base.depwarn(
+            "`Product(v)` is deprecated, please use `product_distribution(v)`",
+            :Product,
+        )
         return new{S, T, V}(v)
     end
 end
@@ -43,26 +48,9 @@ insupport(d::Product, x::AbstractVector) = all(insupport.(d.v, x))
 minimum(d::Product) = map(minimum, d.v)
 maximum(d::Product) = map(maximum, d.v)
 
-"""
-    product_distribution(dists::AbstractVector{<:UnivariateDistribution})
-
-Creates a multivariate product distribution `P` from a vector of univariate distributions.
-Fallback is the `Product constructor`, but specialized methods can be defined
-for distributions with a special multivariate product.
-"""
-function product_distribution(dists::AbstractVector{<:UnivariateDistribution})
-    return Product(dists)
-end
-
-"""
-    product_distribution(dists::AbstractVector{<:Normal})
-
-Computes the multivariate Normal distribution obtained by stacking the univariate
-normal distributions. The result is a multivariate Gaussian with a diagonal
-covariance matrix.
-"""
-function product_distribution(dists::AbstractVector{<:Normal})
-    µ = mean.(dists)
-    σ2 = var.(dists)
-    return MvNormal(µ, Diagonal(σ2))
-end
+# TODO: remove deprecation when `Product` is removed
+# it will return a `ProductDistribution` then which is already the default for
+# higher-dimensional arrays and distributions
+Base.@deprecate product_distribution(
+    dists::AbstractVector{<:UnivariateDistribution}
+) Product(dists)

src/multivariates.jl

@@ -116,7 +116,7 @@ for fname in ["dirichlet.jl",
               "mvnormalcanon.jl",
               "mvlognormal.jl",
               "mvtdist.jl",
-              "product.jl",
+              "product.jl", # deprecated
               "vonmisesfisher.jl"]
     include(joinpath("multivariate", fname))
 end