From c7b9218620607f89a9f96c2035a995cb726984e2 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Thu, 4 Jan 2024 23:49:26 +0100
Subject: [PATCH 1/2] Replace broadcasting over distributions with broadcasting
 with partially applied functions

---
 src/deprecates.jl                             |   8 +-
 src/mixtures/mixturemodel.jl                  |   8 +-
 src/multivariate/jointorderstatistics.jl      |   2 +-
 src/multivariate/mvnormal copy.jl             | 549 ++++++++++++++++++
 src/qq.jl                                     |   4 +-
 src/univariate/continuous/uniform.jl          |   2 +-
 src/univariate/discrete/betabinomial.jl       |  13 +-
 src/univariate/discrete/hypergeometric.jl     |   2 +-
 .../discrete/noncentralhypergeometric.jl      |   9 +-
 test/censored.jl                              |   8 +-
 test/matrixvariates.jl                        |   2 +-
 test/mixture.jl                               |   2 +-
 test/multivariate/product.jl                  |   3 +-
 test/product.jl                               |   2 +-
 test/testutils.jl                             |  46 +-
 test/univariate/continuous/johnsonsu.jl       |   8 +-
 test/univariate/continuous/logitnormal.jl     |   2 +-
 test/univariate/continuous/rician.jl          |  10 +-
 test/univariate/continuous/skewnormal.jl      |   2 +-
 test/univariate/discrete/binomial.jl          |   4 +-
 test/univariate/discrete/categorical.jl       |   4 +-
 test/univariate/discrete/poissonbinomial.jl   |   6 +-
 test/univariate/discrete/soliton.jl           |   4 +-
 test/univariate/locationscale.jl              |   4 +-
 test/univariates.jl                           |   4 +-
 25 files changed, 631 insertions(+), 77 deletions(-)
 create mode 100644 src/multivariate/mvnormal copy.jl

diff --git a/src/deprecates.jl b/src/deprecates.jl
index cacaa12b59..3cfc20a5c6 100644
--- a/src/deprecates.jl
+++ b/src/deprecates.jl
@@ -34,13 +34,13 @@ for fun in [:pdf, :logpdf,
     fun! = Symbol(fun, '!')
 
     @eval begin
-        @deprecate ($_fun!)(r::AbstractArray{<:Real}, d::UnivariateDistribution, X::AbstractArray{<:Real}) r .= ($fun).(d, X) false
-        @deprecate ($fun!)(r::AbstractArray{<:Real}, d::UnivariateDistribution, X::AbstractArray{<:Real}) r .= ($fun).(d, X) false
-        @deprecate ($fun)(d::UnivariateDistribution, X::AbstractArray{<:Real}) ($fun).(d, X)
+        @deprecate ($_fun!)(r::AbstractArray{<:Real}, d::UnivariateDistribution, X::AbstractArray{<:Real}) r .= Base.Fix1($fun, d).(X) false
+        @deprecate ($fun!)(r::AbstractArray{<:Real}, d::UnivariateDistribution, X::AbstractArray{<:Real}) r .= Base.Fix1($fun, d).(X) false
+        @deprecate ($fun)(d::UnivariateDistribution, X::AbstractArray{<:Real}) map(Base.Fix1($fun, d), X)
     end
 end
 
-@deprecate pdf(d::DiscreteUnivariateDistribution) pdf.(Ref(d), support(d))
+@deprecate pdf(d::DiscreteUnivariateDistribution) map(Base.Fix1(pdf, d), support(d))
 
 # Wishart constructors
 @deprecate Wishart(df::Real, S::AbstractPDMat, warn::Bool) Wishart(df, S)
diff --git a/src/mixtures/mixturemodel.jl b/src/mixtures/mixturemodel.jl
index fc05012d0d..c3d3b1f919 100644
--- a/src/mixtures/mixturemodel.jl
+++ b/src/mixtures/mixturemodel.jl
@@ -297,7 +297,7 @@ function _mixpdf!(r::AbstractArray, d::AbstractMixtureModel, x)
         pi = p[i]
         if pi > 0.0
             if d isa UnivariateMixture
-                t .= pdf.(component(d, i), x)
+                t .= Base.Fix1(pdf, component(d, i)).(x)
             else
                 pdf!(t, component(d, i), x)
             end
@@ -326,7 +326,7 @@ function _mixlogpdf!(r::AbstractArray, d::AbstractMixtureModel, x)
             lp_i = view(Lp, :, i)
             # compute logpdf in batch and store
             if d isa UnivariateMixture
-                lp_i .= logpdf.(component(d, i), x)
+                lp_i .= Base.Fix1(logpdf, component(d, i)).(x)
             else
                 logpdf!(lp_i, component(d, i), x)
             end
@@ -398,7 +398,7 @@ function _cwise_pdf!(r::AbstractMatrix, d::AbstractMixtureModel, X)
     size(r) == (n, K) || error("The size of r is incorrect.")
     for i = 1:K
         if d isa UnivariateMixture
-            view(r,:,i) .= pdf.(Ref(component(d, i)), X)
+            view(r,:,i) .= Base.Fix1(pdf, component(d, i)).(X)
         else
             pdf!(view(r,:,i),component(d, i), X)
         end
@@ -412,7 +412,7 @@ function _cwise_logpdf!(r::AbstractMatrix, d::AbstractMixtureModel, X)
     size(r) == (n, K) || error("The size of r is incorrect.")
     for i = 1:K
         if d isa UnivariateMixture
-            view(r,:,i) .= logpdf.(Ref(component(d, i)), X)
+            view(r,:,i) .= Base.Fix1(logpdf, component(d, i)).(X)
         else
             logpdf!(view(r,:,i), component(d, i), X)
         end
diff --git a/src/multivariate/jointorderstatistics.jl b/src/multivariate/jointorderstatistics.jl
index 71c3f93bf3..1fbed0d1b6 100644
--- a/src/multivariate/jointorderstatistics.jl
+++ b/src/multivariate/jointorderstatistics.jl
@@ -162,7 +162,7 @@ function _rand!(rng::AbstractRNG, d::JointOrderStatistics, x::AbstractVector{<:R
         else
             s += randexp(rng, T)
         end
-        x .= quantile.(d.dist, x ./ s)
+        x .= Base.Fix1(quantile, d.dist).(x ./ s)
     end
     return x
 end
diff --git a/src/multivariate/mvnormal copy.jl b/src/multivariate/mvnormal copy.jl
new file mode 100644
index 0000000000..3f0db60f1d
--- /dev/null
+++ b/src/multivariate/mvnormal copy.jl	
@@ -0,0 +1,549 @@
+# Multivariate Normal distribution
+
+
+###########################################################
+#
+#   Abstract base class for multivariate normal
+#
+#   Each subtype should provide the following methods:
+#
+#   - length(d):        vector dimension
+#   - mean(d):          the mean vector (in full form)
+#   - cov(d):           the covariance matrix (in full form)
+#   - invcov(d):        inverse of covariance
+#   - logdetcov(d):     log-determinant of covariance
+#   - sqmahal(d, x):        Squared Mahalanobis distance to center
+#   - sqmahal!(r, d, x):    Squared Mahalanobis distances
+#   - gradlogpdf(d, x):     Gradient of logpdf w.r.t. x
+#   - _rand!(d, x):         Sample random vector(s)
+#
+#   Other generic functions will be implemented on top
+#   of these core functions.
+#
+###########################################################
+
+"""
+
+The [Multivariate normal distribution](http://en.wikipedia.org/wiki/Multivariate_normal_distribution)
+is a multidimensional generalization of the *normal distribution*. The probability density function of
+a d-dimensional multivariate normal distribution with mean vector ``\\boldsymbol{\\mu}`` and
+covariance matrix ``\\boldsymbol{\\Sigma}`` is:
+
+```math
+f(\\mathbf{x}; \\boldsymbol{\\mu}, \\boldsymbol{\\Sigma}) = \\frac{1}{(2 \\pi)^{d/2} |\\boldsymbol{\\Sigma}|^{1/2}}
+\\exp \\left( - \\frac{1}{2} (\\mathbf{x} - \\boldsymbol{\\mu})^T \\Sigma^{-1} (\\mathbf{x} - \\boldsymbol{\\mu}) \\right)
+```
+
+We realize that the mean vector and the covariance often have special forms in practice,
+which can be exploited to simplify the computation. For example, the mean vector is sometimes
+just a zero vector, while the covariance matrix can be a diagonal matrix or even in the form
+of ``\\sigma^2 \\mathbf{I}``. To take advantage of such special cases, we introduce a parametric
+type `MvNormal`, defined as below, which allows users to specify the special structure of
+the mean and covariance.
+
+```julia
+struct MvNormal{T<:Real,Cov<:AbstractPDMat,Mean<:AbstractVector} <: AbstractMvNormal
+    μ::Mean
+    Σ::Cov
+end
+```
+
+Here, the mean vector can be an instance of any `AbstractVector`. The covariance can be
+of any subtype of `AbstractPDMat`. Particularly, one can use `PDMat` for full covariance,
+`PDiagMat` for diagonal covariance, and `ScalMat` for the isotropic covariance -- those
+in the form of ``\\sigma^2 \\mathbf{I}``. (See the Julia package
+[PDMats](https://github.com/JuliaStats/PDMats.jl/) for details).
+
+We also define a set of aliases for the types using different combinations of mean vectors and covariance:
+
+```julia
+const IsoNormal  = MvNormal{Float64, ScalMat{Float64},                  Vector{Float64}}
+const DiagNormal = MvNormal{Float64, PDiagMat{Float64,Vector{Float64}}, Vector{Float64}}
+const FullNormal = MvNormal{Float64, PDMat{Float64,Matrix{Float64}},    Vector{Float64}}
+
+const ZeroMeanIsoNormal{Axes}  = MvNormal{Float64, ScalMat{Float64},                  Zeros{Float64,1,Axes}}
+const ZeroMeanDiagNormal{Axes} = MvNormal{Float64, PDiagMat{Float64,Vector{Float64}}, Zeros{Float64,1,Axes}}
+const ZeroMeanFullNormal{Axes} = MvNormal{Float64, PDMat{Float64,Matrix{Float64}},    Zeros{Float64,1,Axes}}
+```
+
+Multivariate normal distributions support affine transformations:
+```julia
+d = MvNormal(μ, Σ)
+c + B * d    # == MvNormal(B * μ + c, B * Σ * B')
+dot(b, d)    # == Normal(dot(b, μ), b' * Σ * b)
+```
+"""
+abstract type AbstractMvNormal <: ContinuousMultivariateDistribution end
+
+### Generic methods (for all AbstractMvNormal subtypes)
+
+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))
+mode(d::AbstractMvNormal) = mean(d)
+modes(d::AbstractMvNormal) = [mean(d)]
+
+"""
+    rand(::AbstractRNG, ::Distributions.AbstractMvNormal)
+
+Sample a random vector from the provided multi-variate normal distribution.
+"""
+rand(::AbstractRNG, ::Distributions.AbstractMvNormal)
+
+function entropy(d::AbstractMvNormal)
+    ldcd = logdetcov(d)
+    return (length(d) * (oftype(ldcd, log2π) + 1) + ldcd) / 2
+end
+
+function mvnormal_c0(d::AbstractMvNormal)
+    ldcd = logdetcov(d)
+    return - (length(d) * oftype(ldcd, log2π) + ldcd) / 2
+end
+
+function kldivergence(p::AbstractMvNormal, q::AbstractMvNormal)
+    # This is the generic implementation for AbstractMvNormal, you might need to specialize for your type
+    length(p) == length(q) ||
+        throw(DimensionMismatch("Distributions p and q have different dimensions $(length(p)) and $(length(q))"))
+    # logdetcov is used separately from _cov for any potential optimization done there
+    return (tr(_cov(q) \ _cov(p)) + sqmahal(q, mean(p)) - length(p) + logdetcov(q) - logdetcov(p)) / 2
+end
+
+# This is a workaround to take advantage of the PDMats objects for MvNormal and avoid copies as Matrix
+# TODO: Remove this once `cov(::MvNormal)` returns the PDMats object
+_cov(d::AbstractMvNormal) = cov(d)
+
+"""
+    invcov(d::AbstractMvNormal)
+
+Return the inversed covariance matrix of d.
+"""
+invcov(d::AbstractMvNormal)
+
+"""
+    logdetcov(d::AbstractMvNormal)
+
+Return the log-determinant value of the covariance matrix.
+"""
+logdetcov(d::AbstractMvNormal)
+
+"""
+    sqmahal(d::AbstractMvNormal, x::AbstractVecOrMat)
+
+Return the squared Mahalanobis distance from `x` to the mean of `d` w.r.t. the covariance.
+When x is a vector, it returns a scalar value. When x is a matrix, it returns a vector of length size(x,2).
+"""
+sqmahal(d::AbstractMvNormal, x::AbstractVecOrMat)
+
+"""
+    sqmahal!(r::AbstractVector, d::AbstractMvNormal, x::AbstractMatrix)
+
+Compute [`sqmahal(d, x)`](@ref) and write the results to the pre-allocated vector `r`.
+"""
+sqmahal!(r::AbstractVector, d::AbstractMvNormal, x::AbstractMatrix)
+
+@inline function logpdf(d::AbstractMvNormal, x::AbstractArray{<:Real,N}) where {N}
+    @boundscheck begin
+        size(x, 1) == length(d) || throw(DimensionMismatch("inconsistent array dimensions"))
+    end
+
+    # Compute constants (with respect to `x`) only once
+    c0 = mvnormal_c0(d) 
+
+    if N == 1
+        # Just numbers in this case
+        return c0 - sqmahal(d, x) / 2
+    else
+        # Compute squared Mahalanobis distance
+        xmat = N == 2 ? x : reshape(x, size(x, 1), :)
+        sqm = sqmahal(d, xmat)
+        
+        # Optimization with reduced allocations for Array
+        res = if sqm isa Array
+            @inbounds for i in eachindex(sqm)
+                sqm[i] = c0 - sqm[i] / 2
+            end
+            sqm
+        else
+            c0 .- sqm ./ 2
+        end
+
+        return N == 2 ? res : reshape(res, ntuple(i -> size(x, i + 1), Val(N - 1)))
+    end
+end
+
+function _logpdf!(r::AbstractArray{<:Real}, d::AbstractMvNormal, x::AbstractArray{<:Real,N}) where {N}
+    dim = length(d)
+    xmat = N == 2 ? x : reshape(x, dim, :)
+    sqmahal!(vec(r), d, xmat)
+    c0 = mvnormal_c0(d)
+    for i in eachindex(r)
+        @inbounds r[i] = c0 - r[i]/2
+    end
+    return r
+end
+
+###########################################################
+#
+#   MvNormal
+#
+#   Multivariate normal distribution with mean parameters
+#
+###########################################################
+"""
+    MvNormal
+
+Generally, users don't have to worry about these internal details.
+
+We provide a common constructor `MvNormal`, which will construct a distribution of
+appropriate type depending on the input arguments.
+"""
+struct MvNormal{T<:Real,Cov<:AbstractPDMat,Mean<:AbstractVector} <: AbstractMvNormal
+    μ::Mean
+    Σ::Cov
+end
+
+const MultivariateNormal = MvNormal  # for the purpose of backward compatibility
+
+const IsoNormal  = MvNormal{Float64,ScalMat{Float64},Vector{Float64}}
+const DiagNormal = MvNormal{Float64,PDiagMat{Float64,Vector{Float64}},Vector{Float64}}
+const FullNormal = MvNormal{Float64,PDMat{Float64,Matrix{Float64}},Vector{Float64}}
+
+const ZeroMeanIsoNormal{Axes}  = MvNormal{Float64,ScalMat{Float64},Zeros{Float64,1,Axes}}
+const ZeroMeanDiagNormal{Axes} = MvNormal{Float64,PDiagMat{Float64,Vector{Float64}},Zeros{Float64,1,Axes}}
+const ZeroMeanFullNormal{Axes} = MvNormal{Float64,PDMat{Float64,Matrix{Float64}},Zeros{Float64,1,Axes}}
+
+### Construction
+function MvNormal(μ::AbstractVector{T}, Σ::AbstractPDMat{T}) where {T<:Real}
+    size(Σ, 1) == length(μ) || throw(DimensionMismatch("The dimensions of mu and Sigma are inconsistent."))
+    MvNormal{T,typeof(Σ), typeof(μ)}(μ, Σ)
+end
+
+function MvNormal(μ::AbstractVector{<:Real}, Σ::AbstractPDMat{<:Real})
+    R = Base.promote_eltype(μ, Σ)
+    MvNormal(convert(AbstractArray{R}, μ), convert(AbstractArray{R}, Σ))
+end
+
+# constructor with general covariance matrix
+"""
+    MvNormal(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real})
+
+Construct a multivariate normal distribution with mean `μ` and covariance matrix `Σ`.
+"""
+MvNormal(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real}) = MvNormal(μ, PDMat(Σ))
+MvNormal(μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real}) = MvNormal(μ, PDiagMat(Σ.diag))
+MvNormal(μ::AbstractVector{<:Real}, Σ::Union{Symmetric{<:Real,<:Diagonal{<:Real}},Hermitian{<:Real,<:Diagonal{<:Real}}}) = MvNormal(μ, PDiagMat(Σ.data.diag))
+MvNormal(μ::AbstractVector{<:Real}, Σ::UniformScaling{<:Real}) =
+    MvNormal(μ, ScalMat(length(μ), Σ.λ))
+function MvNormal(
+    μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real,<:FillArrays.AbstractFill{<:Real,1}}
+)
+    return MvNormal(μ, ScalMat(size(Σ, 1), FillArrays.getindex_value(Σ.diag)))
+end
+
+# constructor without mean vector
+"""
+    MvNormal(Σ::AbstractMatrix{<:Real})
+
+Construct a multivariate normal distribution with zero mean and covariance matrix `Σ`.
+"""
+MvNormal(Σ::AbstractMatrix{<:Real}) = MvNormal(Zeros{eltype(Σ)}(size(Σ, 1)), Σ)
+
+# deprecated constructors with standard deviations
+Base.@deprecate MvNormal(μ::AbstractVector{<:Real}, σ::AbstractVector{<:Real}) MvNormal(μ, LinearAlgebra.Diagonal(map(abs2, σ)))
+Base.@deprecate MvNormal(μ::AbstractVector{<:Real}, σ::Real) MvNormal(μ, σ^2 * I)
+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.Σ))
+end
+Base.convert(::Type{MvNormal{T}}, d::MvNormal{T}) where {T<:Real} = d
+
+function convert(::Type{MvNormal{T}}, μ::AbstractVector, Σ::AbstractPDMat) where T<:Real
+    MvNormal(convert(AbstractArray{T}, μ), convert(AbstractArray{T}, Σ))
+end
+
+### Show
+
+distrname(d::IsoNormal)  = "IsoNormal"    # Note: IsoNormal, etc are just alias names.
+distrname(d::DiagNormal) = "DiagNormal"
+distrname(d::FullNormal) = "FullNormal"
+
+distrname(d::ZeroMeanIsoNormal) = "ZeroMeanIsoNormal"
+distrname(d::ZeroMeanDiagNormal) = "ZeroMeanDiagNormal"
+distrname(d::ZeroMeanFullNormal) = "ZeroMeanFullNormal"
+
+Base.show(io::IO, d::MvNormal) =
+    show_multline(io, d, [(:dim, length(d)), (:μ, mean(d)), (:Σ, cov(d))])
+
+### Basic statistics
+
+length(d::MvNormal) = length(d.μ)
+mean(d::MvNormal) = d.μ
+params(d::MvNormal) = (d.μ, d.Σ)
+@inline partype(d::MvNormal{T}) where {T<:Real} = T
+
+var(d::MvNormal) = diag(d.Σ)
+cov(d::MvNormal) = Matrix(d.Σ)
+_cov(d::MvNormal) = d.Σ
+
+invcov(d::MvNormal) = Matrix(inv(d.Σ))
+logdetcov(d::MvNormal) = logdet(d.Σ)
+
+### Evaluation
+
+sqmahal(d::MvNormal, x::AbstractVector) = invquad(d.Σ, x - d.μ)
+sqmahal(d::MvNormal, x::AbstractMatrix) = invquad(d.Σ, x .- d.μ)
+sqmahal!(r::AbstractVector, d::MvNormal, x::AbstractMatrix) = invquad!(r, d.Σ, x .- d.μ)
+
+gradlogpdf(d::MvNormal, x::AbstractVector{<:Real}) = -(d.Σ \ (x .- d.μ))
+
+# Sampling (for GenericMvNormal)
+
+function _rand!(rng::AbstractRNG, d::MvNormal, x::VecOrMat)
+    unwhiten!(d.Σ, randn!(rng, x))
+    x .+= d.μ
+    return x
+end
+
+# Workaround: randn! only works for Array, but not generally for AbstractArray
+function _rand!(rng::AbstractRNG, d::MvNormal, x::AbstractVector)
+    for i in eachindex(x)
+        @inbounds x[i] = randn(rng, eltype(x))
+    end
+    unwhiten!(d.Σ, x)
+    x .+= d.μ
+    return x
+end
+
+### Affine transformations
+
+Base.:+(d::MvNormal, c::AbstractVector) = MvNormal(d.μ + c, d.Σ)
+Base.:+(c::AbstractVector, d::MvNormal) = d + c
+Base.:-(d::MvNormal, c::AbstractVector) = MvNormal(d.μ - c, d.Σ)
+
+Base.:*(B::AbstractMatrix, d::MvNormal) = MvNormal(B * d.μ, X_A_Xt(d.Σ, B))
+
+dot(b::AbstractVector, d::MvNormal) = Normal(dot(d.μ, b), √quad(d.Σ, b))
+
+dot(d::MvNormal, b::AbstractVector) = dot(b, d)
+
+###########################################################
+#
+#   Estimation of MvNormal
+#
+###########################################################
+
+### Estimation with known covariance
+
+struct MvNormalKnownCov{Cov<:AbstractPDMat}
+    Σ::Cov
+end
+
+MvNormalKnownCov(d::Int, σ::Real) = MvNormalKnownCov(ScalMat(d, abs2(Float64(σ))))
+MvNormalKnownCov(σ::Vector{Float64}) = MvNormalKnownCov(PDiagMat(abs2.(σ)))
+MvNormalKnownCov(Σ::Matrix{Float64}) = MvNormalKnownCov(PDMat(Σ))
+
+length(g::MvNormalKnownCov) = size(g.Σ, 1)
+
+struct MvNormalKnownCovStats{Cov<:AbstractPDMat}
+    invΣ::Cov              # inverse covariance
+    sx::Vector{Float64}    # (weighted) sum of vectors
+    tw::Float64            # sum of weights
+end
+
+function suffstats(g::MvNormalKnownCov{Cov}, x::AbstractMatrix{Float64}) where Cov<:AbstractPDMat
+    size(x,1) == length(g) || throw(DimensionMismatch("Invalid argument dimensions."))
+    invΣ = inv(g.Σ)
+    sx = vec(sum(x, dims=2))
+    tw = Float64(size(x, 2))
+    MvNormalKnownCovStats{Cov}(invΣ, sx, tw)
+end
+
+function suffstats(g::MvNormalKnownCov{Cov}, x::AbstractMatrix{Float64}, w::AbstractVector) where Cov<:AbstractPDMat
+    (size(x,1) == length(g) && size(x,2) == length(w)) ||
+        throw(DimensionMismatch("Inconsistent argument dimensions."))
+    invΣ = inv(g.Σ)
+    sx = x * vec(w)
+    tw = sum(w)
+    MvNormalKnownCovStats{Cov}(invΣ, sx, tw)
+end
+
+## MLE estimation with covariance known
+
+fit_mle(g::MvNormalKnownCov{C}, ss::MvNormalKnownCovStats{C}) where {C<:AbstractPDMat} =
+    MvNormal(ss.sx * inv(ss.tw), g.Σ)
+
+function fit_mle(g::MvNormalKnownCov, x::AbstractMatrix{Float64})
+    d = length(g)
+    size(x,1) == d || throw(DimensionMismatch("Invalid argument dimensions."))
+    μ = lmul!(inv(size(x,2)), vec(sum(x,dims=2)))
+    MvNormal(μ, g.Σ)
+end
+
+function fit_mle(g::MvNormalKnownCov, x::AbstractMatrix{Float64}, w::AbstractVector)
+    d = length(g)
+    (size(x,1) == d && size(x,2) == length(w)) ||
+        throw(DimensionMismatch("Inconsistent argument dimensions."))
+    μ = BLAS.gemv('N', inv(sum(w)), x, vec(w))
+    MvNormal(μ, g.Σ)
+end
+
+
+### Estimation (both mean and cov unknown)
+
+struct MvNormalStats <: SufficientStats
+    s::Vector{Float64}  # (weighted) sum of x
+    m::Vector{Float64}  # (weighted) mean of x
+    s2::Matrix{Float64} # (weighted) sum of (x-μ) * (x-μ)'
+    tw::Float64         # total sample weight
+end
+
+function suffstats(D::Type{MvNormal}, x::AbstractMatrix{Float64})
+    d = size(x, 1)
+    n = size(x, 2)
+    s = vec(sum(x, dims=2))
+    m = s * inv(n)
+    z = x .- m
+    s2 = z * z'
+    MvNormalStats(s, m, s2, Float64(n))
+end
+
+function suffstats(D::Type{MvNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
+    d = size(x, 1)
+    n = size(x, 2)
+    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions."))
+
+    tw = sum(w)
+    s = x * vec(w)
+    m = s * inv(tw)
+    z = similar(x)
+    for j = 1:n
+        xj = view(x,:,j)
+        zj = view(z,:,j)
+        swj = sqrt(w[j])
+        for i = 1:d
+            @inbounds zj[i] = swj * (xj[i] - m[i])
+        end
+    end
+    s2 = z * z'
+    MvNormalStats(s, m, s2, tw)
+end
+
+
+# Maximum Likelihood Estimation
+#
+# Specialized algorithms are respectively implemented for
+# each kind of covariance
+#
+
+fit_mle(D::Type{MvNormal}, ss::MvNormalStats) = fit_mle(FullNormal, ss)
+fit_mle(D::Type{MvNormal}, x::AbstractMatrix{Float64}) = fit_mle(FullNormal, x)
+fit_mle(D::Type{MvNormal}, x::AbstractMatrix{Float64}, w::AbstractArray{Float64}) = fit_mle(FullNormal, x, w)
+
+fit_mle(D::Type{FullNormal}, ss::MvNormalStats) = MvNormal(ss.m, ss.s2 * inv(ss.tw))
+
+function fit_mle(D::Type{FullNormal}, x::AbstractMatrix{Float64})
+    n = size(x, 2)
+    mu = vec(mean(x, dims=2))
+    z = x .- mu
+    C = BLAS.syrk('U', 'N', inv(n), z)
+    LinearAlgebra.copytri!(C, 'U')
+    MvNormal(mu, PDMat(C))
+end
+
+function fit_mle(D::Type{FullNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
+    m = size(x, 1)
+    n = size(x, 2)
+    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions"))
+
+    inv_sw = inv(sum(w))
+    mu = BLAS.gemv('N', inv_sw, x, w)
+
+    z = Matrix{Float64}(undef, m, n)
+    for j = 1:n
+        cj = sqrt(w[j])
+        for i = 1:m
+            @inbounds z[i,j] = (x[i,j] - mu[i]) * cj
+        end
+    end
+    C = BLAS.syrk('U', 'N', inv_sw, z)
+    LinearAlgebra.copytri!(C, 'U')
+    MvNormal(mu, PDMat(C))
+end
+
+function fit_mle(D::Type{DiagNormal}, x::AbstractMatrix{Float64})
+    m = size(x, 1)
+    n = size(x, 2)
+
+    mu = vec(mean(x, dims=2))
+    va = zeros(Float64, m)
+    for j = 1:n
+        for i = 1:m
+            @inbounds va[i] += abs2(x[i,j] - mu[i])
+        end
+    end
+    lmul!(inv(n), va)
+    MvNormal(mu, PDiagMat(va))
+end
+
+function fit_mle(D::Type{DiagNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
+    m = size(x, 1)
+    n = size(x, 2)
+    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions"))
+
+    inv_sw = inv(sum(w))
+    mu = BLAS.gemv('N', inv_sw, x, w)
+
+    va = zeros(Float64, m)
+    for j = 1:n
+        @inbounds wj = w[j]
+        for i = 1:m
+            @inbounds va[i] += abs2(x[i,j] - mu[i]) * wj
+        end
+    end
+    lmul!(inv_sw, va)
+    MvNormal(mu, PDiagMat(va))
+end
+
+function fit_mle(D::Type{IsoNormal}, x::AbstractMatrix{Float64})
+    m = size(x, 1)
+    n = size(x, 2)
+
+    mu = vec(mean(x, dims=2))
+    va = 0.
+    for j = 1:n
+        va_j = 0.
+        for i = 1:m
+            @inbounds va_j += abs2(x[i,j] - mu[i])
+        end
+        va += va_j
+    end
+    MvNormal(mu, ScalMat(m, va / (m * n)))
+end
+
+function fit_mle(D::Type{IsoNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
+    m = size(x, 1)
+    n = size(x, 2)
+    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions"))
+
+    sw = sum(w)
+    inv_sw = 1.0 / sw
+    mu = BLAS.gemv('N', inv_sw, x, w)
+
+    va = 0.
+    for j = 1:n
+        @inbounds wj = w[j]
+        va_j = 0.
+        for i = 1:m
+            @inbounds va_j += abs2(x[i,j] - mu[i]) * wj
+        end
+        va += va_j
+    end
+    MvNormal(mu, ScalMat(m, va / (m * sw)))
+end
diff --git a/src/qq.jl b/src/qq.jl
index fe80e31a68..6431ed2bf9 100644
--- a/src/qq.jl
+++ b/src/qq.jl
@@ -35,14 +35,14 @@ function qqbuild(x::AbstractVector, d::UnivariateDistribution)
     n = length(x)
     grid = ppoints(n)
     qx = quantile(x, grid)
-    qd = quantile.(Ref(d), grid)
+    qd = map(Base.Fix1(quantile, d), grid)
     return QQPair(qx, qd)
 end
 
 function qqbuild(d::UnivariateDistribution, x::AbstractVector)
     n = length(x)
     grid = ppoints(n)
-    qd = quantile.(Ref(d), grid)
+    qd = map(Base.Fix1(quantile, d), grid)
     qx = quantile(x, grid)
     return QQPair(qd, qx)
 end
diff --git a/src/univariate/continuous/uniform.jl b/src/univariate/continuous/uniform.jl
index 1f535159d0..d4beca7d79 100644
--- a/src/univariate/continuous/uniform.jl
+++ b/src/univariate/continuous/uniform.jl
@@ -154,7 +154,7 @@ Base.:*(c::Real, d::Uniform) = Uniform(minmax(c * d.a, c * d.b)...)
 rand(rng::AbstractRNG, d::Uniform) = d.a + (d.b - d.a) * rand(rng)
 
 _rand!(rng::AbstractRNG, d::Uniform, A::AbstractArray{<:Real}) =
-    A .= quantile.(d, rand!(rng, A))
+    A .= Base.Fix1(quantile, d).(rand!(rng, A))
 
 
 #### Fitting
diff --git a/src/univariate/discrete/betabinomial.jl b/src/univariate/discrete/betabinomial.jl
index 618b3f83ef..7f2a0ac0dc 100644
--- a/src/univariate/discrete/betabinomial.jl
+++ b/src/univariate/discrete/betabinomial.jl
@@ -103,12 +103,15 @@ for f in (:ccdf, :logcdf, :logccdf)
     end
 end
 
-entropy(d::BetaBinomial) = entropy(Categorical(pdf.(Ref(d),support(d))))
-median(d::BetaBinomial) = median(Categorical(pdf.(Ref(d),support(d)))) - 1
-mode(d::BetaBinomial) = argmax(pdf.(Ref(d),support(d))) - 1
-modes(d::BetaBinomial) = modes(Categorical(pdf.(Ref(d),support(d)))) .- 1
+# Shifted categorical distribution corresponding to `BetaBinomial`
+_categorical(d::BetaBinomial) = Categorical(map(Base.Fix1(pdf, d), support(d)))
 
-quantile(d::BetaBinomial, p::Float64) = quantile(Categorical(pdf.(Ref(d), support(d))), p) - 1
+entropy(d::BetaBinomial) = entropy(_categorical(d))
+median(d::BetaBinomial) = median(_categorical(d)) - 1
+mode(d::BetaBinomial) = mode(_categorical(d)) - 1
+modes(d::BetaBinomial) = modes(_categorical(d)) .- 1
+
+quantile(d::BetaBinomial, p::Float64) = quantile(_categorical(d), p) - 1
 
 #### Sampling
 
diff --git a/src/univariate/discrete/hypergeometric.jl b/src/univariate/discrete/hypergeometric.jl
index 9f72c2fea4..de575abdcb 100644
--- a/src/univariate/discrete/hypergeometric.jl
+++ b/src/univariate/discrete/hypergeometric.jl
@@ -75,7 +75,7 @@ function kurtosis(d::Hypergeometric)
     a/b
 end
 
-entropy(d::Hypergeometric) = entropy(pdf.(Ref(d), support(d)))
+entropy(d::Hypergeometric) = entropy(map(Base.Fix1(pdf, d), support(d)))
 
 ### Evaluation & Sampling
 
diff --git a/src/univariate/discrete/noncentralhypergeometric.jl b/src/univariate/discrete/noncentralhypergeometric.jl
index 4c33d26215..52853cfc08 100644
--- a/src/univariate/discrete/noncentralhypergeometric.jl
+++ b/src/univariate/discrete/noncentralhypergeometric.jl
@@ -256,9 +256,12 @@ end
 Base.convert(::Type{WalleniusNoncentralHypergeometric{T}}, d::WalleniusNoncentralHypergeometric{T}) where {T<:Real} = d
 
 # Properties
-mean(d::WalleniusNoncentralHypergeometric) = sum(support(d) .* pdf.(Ref(d), support(d)))
-var(d::WalleniusNoncentralHypergeometric)  = sum((support(d) .- mean(d)).^2 .* pdf.(Ref(d), support(d)))
-mode(d::WalleniusNoncentralHypergeometric) = support(d)[argmax(pdf.(Ref(d), support(d)))]
+function _discretenonparametric(d::WalleniusNoncentralHypergeometric)
+    return DiscreteNonParametric(support(d), map(Base.Fix1(pdf, d), support(d)))
+end
+mean(d::WalleniusNoncentralHypergeometric) = mean(_discretenonparametric(d))
+var(d::WalleniusNoncentralHypergeometric)  = var(_discretenonparametric(d))
+mode(d::WalleniusNoncentralHypergeometric) = mode(_discretenonparametric(d))
 
 entropy(d::WalleniusNoncentralHypergeometric) = 1
 
diff --git a/test/censored.jl b/test/censored.jl
index eaad72cdc2..27f30cadfc 100644
--- a/test/censored.jl
+++ b/test/censored.jl
@@ -207,7 +207,7 @@ end
             end
             @test @inferred(median(d)) ≈ clamp(median(d0), l, u)
             @inferred quantile(d, 0.5)
-            @test quantile.(d, 0:0.01:1) ≈ clamp.(quantile.(d0, 0:0.01:1), l, u)
+            @test Base.Fix1(quantile, d).(0:0.01:1) ≈ clamp.(Base.Fix1(quantile, d0).(0:0.01:1), l, u)
             # special-case pdf/logpdf/loglikelihood since when replacing Dirac(μ) with
             # Normal(μ, 0), they are infinite
             if lower === nothing || !isfinite(lower)
@@ -253,7 +253,7 @@ end
                 @test f(d) ≈ f(dmix)
             end
             @test median(d) ≈ clamp(median(d0), l, u)
-            @test quantile.(d, 0:0.01:1) ≈ clamp.(quantile.(d0, 0:0.01:1), l, u)
+            @test Base.Fix1(quantile, d).(0:0.01:1) ≈ clamp.(Base.Fix1(quantile, d0).(0:0.01:1), l, u)
             # special-case pdf/logpdf/loglikelihood since when replacing Dirac(μ) with
             # Normal(μ, 0), they are infinite
             if lower === nothing
@@ -311,7 +311,7 @@ end
             end
             @test @inferred(median(d)) ≈ clamp(median(d0), l, u)
             @inferred quantile(d, 0.5)
-            @test quantile.(d, 0:0.01:1) ≈ clamp.(quantile.(d0, 0:0.01:1), l, u)
+            @test Base.Fix1(quantile, d).(0:0.01:1) ≈ clamp.(Base.Fix1(quantile, d0).(0:0.01:1), l, u)
             # rand
             x = rand(d, 10_000)
             @test all(x -> insupport(d, x), x)
@@ -346,7 +346,7 @@ end
                 @test f(d, 5) ≈ f(dmix, 5)
             end
             @test median(d) ≈ clamp(median(d0), l, u)
-            @test quantile.(d, 0:0.01:0.99) ≈ clamp.(quantile.(d0, 0:0.01:0.99), l, u)
+            @test Base.Fix1(quantile, d).(0:0.01:0.99) ≈ clamp.(Base.Fix1(quantile, d0).(0:0.01:0.99), l, u)
             x = rand(d, 100)
             @test loglikelihood(d, x) ≈ loglikelihood(dmix, x)
             # rand
diff --git a/test/matrixvariates.jl b/test/matrixvariates.jl
index 3d0dddfff6..18e984a470 100644
--- a/test/matrixvariates.jl
+++ b/test/matrixvariates.jl
@@ -117,7 +117,7 @@ function test_convert(d::MatrixDistribution)
     @test d == deepcopy(d)
     for elty in (Float32, Float64, BigFloat)
         del1 = convert(distname{elty}, d)
-        del2 = convert(distname{elty}, getfield.(Ref(d), fieldnames(typeof(d)))...)
+        del2 = convert(distname{elty}, (Base.Fix1(getfield, d)).(fieldnames(typeof(d)))...)
         @test del1 isa distname{elty}
         @test del2 isa distname{elty}
         @test partype(del1) == elty
diff --git a/test/mixture.jl b/test/mixture.jl
index 9b84f69326..0b25a2346a 100644
--- a/test/mixture.jl
+++ b/test/mixture.jl
@@ -40,7 +40,7 @@ function test_mixture(g::UnivariateMixture, n::Int, ns::Int,
     for i = 1:n
         @test @inferred(cdf(g, X[i])) ≈ cf[i]
     end
-    @test cdf.(g, X) ≈ cf
+    @test Base.Fix1(cdf, g).(X) ≈ cf
 
     # evaluation
     P0 = zeros(T, n, K)
diff --git a/test/multivariate/product.jl b/test/multivariate/product.jl
index 840eb409f2..c452f96f76 100644
--- a/test/multivariate/product.jl
+++ b/test/multivariate/product.jl
@@ -67,8 +67,7 @@ end
 
     for a in ([0, 1], [-0.5, 0.5])
         # Construct independent distributions and `Product` distribution from these.
-        support = fill(a, N)
-        ds = DiscreteNonParametric.(support, Ref([0.5, 0.5]))
+        ds = [DiscreteNonParametric(copy(a), [0.5, 0.5]) for _ in 1:N]
         x = rand.(ds)
         d_product = product_distribution(ds)
         @test d_product isa Product
diff --git a/test/product.jl b/test/product.jl
index 829e80d1b8..16f12328d3 100644
--- a/test/product.jl
+++ b/test/product.jl
@@ -93,7 +93,7 @@ end
 
     for a in ([0, 1], [-0.5, 0.5])
         # Construct independent distributions and `ProductDistribution` from these.
-        ds1 = DiscreteNonParametric.(fill(a, N), Ref([0.5, 0.5]))
+        ds1 = [DiscreteNonParametric(copy(a), [0.5, 0.5]) for _ in 1:N]
         # Replace with
         # d_product1 = @inferred(product_distribution(ds1))
         # when `Product` is removed
diff --git a/test/testutils.jl b/test/testutils.jl
index ad109e58e3..2859856a73 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -128,7 +128,7 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
     rmin = floor(Int,quantile(distr, 0.00001))::Int
     rmax = floor(Int,quantile(distr, 0.99999))::Int
     m = rmax - rmin + 1  # length of the range
-    p0 = pdf.(Ref(distr), rmin:rmax)  # reference probability masses
+    p0 = map(Base.Fix1(pdf, distr), rmin:rmax)  # reference probability masses
     @assert length(p0) == m
 
     # determine confidence intervals for counts:
@@ -233,7 +233,7 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
     #
     clb = Vector{Int}(undef, nbins)
     cub = Vector{Int}(undef, nbins)
-    cdfs = cdf.(Ref(distr), edges)
+    cdfs = map(Base.Fix1(cdf, distr), edges)
 
     for i = 1:nbins
         pi = cdfs[i+1] - cdfs[i]
@@ -385,19 +385,19 @@ function test_range_evaluation(d::DiscreteUnivariateDistribution)
     rmin = round(Int, islowerbounded(d) ? vmin : quantile(d, 0.001))::Int
     rmax = round(Int, isupperbounded(d) ? vmax : quantile(d, 0.999))::Int
 
-    p0 = pdf.(Ref(d), collect(rmin:rmax))
-    @test pdf.(Ref(d), rmin:rmax) ≈ p0
+    p0 = map(Base.Fix1(pdf, d), collect(rmin:rmax))
+    @test map(Base.Fix1(pdf, d), rmin:rmax) ≈ p0
     if rmin + 2 <= rmax
-        @test pdf.(Ref(d), rmin+1:rmax-1) ≈ p0[2:end-1]
+        @test map(Base.Fix1(pdf, d), rmin+1:rmax-1) ≈ p0[2:end-1]
     end
 
     if isbounded(d)
-        @test pdf.(Ref(d), support(d)) ≈ p0
-        @test pdf.(Ref(d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0)
-        @test pdf.(Ref(d), rmin:rmax+3) ≈ vcat(p0, 0.0, 0.0, 0.0)
-        @test pdf.(Ref(d), rmin-2:rmax+3) ≈ vcat(0.0, 0.0, p0, 0.0, 0.0, 0.0)
+        @test map(Base.Fix1(pdf, d), support(d)) ≈ p0
+        @test map(Base.Fix1(pdf, d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0)
+        @test map(Base.Fix1(pdf, d), rmin:rmax+3) ≈ vcat(p0, 0.0, 0.0, 0.0)
+        @test map(Base.Fix1(pdf, d), rmin-2:rmax+3) ≈ vcat(0.0, 0.0, p0, 0.0, 0.0, 0.0)
     elseif islowerbounded(d)
-        @test pdf.(Ref(d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0)
+        @test map(Base.Fix1(pdf, d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0)
     end
 end
 
@@ -444,13 +444,13 @@ function test_evaluation(d::DiscreteUnivariateDistribution, vs::AbstractVector,
     end
 
     # consistency of scalar-based and vectorized evaluation
-    @test pdf.(Ref(d), vs)  ≈ p
-    @test cdf.(Ref(d), vs)  ≈ c
-    @test ccdf.(Ref(d), vs) ≈ cc
+    @test Base.Fix1(pdf, d).(vs)  ≈ p
+    @test Base.Fix1(cdf, d).(vs)  ≈ c
+    @test Base.Fix1(ccdf, d).(vs) ≈ cc
 
-    @test logpdf.(Ref(d), vs)  ≈ lp
-    @test logcdf.(Ref(d), vs)  ≈ lc
-    @test logccdf.(Ref(d), vs) ≈ lcc
+    @test Base.Fix1(logpdf, d).(vs)  ≈ lp
+    @test Base.Fix1(logcdf, d).(vs)  ≈ lc
+    @test Base.Fix1(logccdf, d).(vs) ≈ lcc
 end
 
 
@@ -511,15 +511,15 @@ function test_evaluation(d::ContinuousUnivariateDistribution, vs::AbstractVector
 
     # consistency of scalar-based and vectorized evaluation
     if !isa(d, StudentizedRange)
-        @test pdf.(Ref(d), vs)  ≈ p
-        @test logpdf.(Ref(d), vs)  ≈ lp
+        @test Base.Fix1(pdf, d).(vs) ≈ p
+        @test Base.Fix1(logpdf, d).(vs) ≈ lp
     end
 
-    @test cdf.(Ref(d), vs)  ≈ c
-    @test ccdf.(Ref(d), vs) ≈ cc
+    @test Base.Fix1(cdf, d).(vs)  ≈ c
+    @test Base.Fix1(ccdf, d).(vs) ≈ cc
 
-    @test logcdf.(Ref(d), vs)  ≈ lc
-    @test logccdf.(Ref(d), vs) ≈ lcc
+    @test Base.Fix1(logcdf, d).(vs)  ≈ lc
+    @test Base.Fix1(logccdf, d).(vs) ≈ lcc
 end
 
 function test_nonfinite(distr::UnivariateDistribution)
@@ -550,7 +550,7 @@ function test_stats(d::DiscreteUnivariateDistribution, vs::AbstractVector)
     # using definition (or an approximation)
 
     vf = Float64[v for v in vs]
-    p = pdf.(Ref(d), vf)
+    p = Base.Fix1(pdf, d).(vf)
     xmean = dot(p, vf)
     xvar = dot(p, abs2.(vf .- xmean))
     xstd = sqrt(xvar)
diff --git a/test/univariate/continuous/johnsonsu.jl b/test/univariate/continuous/johnsonsu.jl
index 5f259167f8..716f1b1df8 100644
--- a/test/univariate/continuous/johnsonsu.jl
+++ b/test/univariate/continuous/johnsonsu.jl
@@ -10,10 +10,10 @@
     @test rand(d1) isa Float64
 
     @test median(d1) == quantile(d1, 0.5)
-    x = quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90])
-    @test all(cdf.(d1, x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90])
-    y = cquantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90])
-    @test all(ccdf.(d1, y) .≈ [0.25, 0.45, 0.60, 0.80, 0.90])
+    x = Base.Fix1(quantile, d1).([0.25, 0.45, 0.60, 0.80, 0.90])
+    @test all(Base.Fix1(cdf, d1).(x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90])
+    y = Base.Fix1(cquantile, d1).([0.25, 0.45, 0.60, 0.80, 0.90])
+    @test all(Base.Fix1(ccdf, d1).(y) .≈ [0.25, 0.45, 0.60, 0.80, 0.90])
 
     @test mean(d1) ≈ 7.581281
     @test var(d1) ≈ 19.1969485
diff --git a/test/univariate/continuous/logitnormal.jl b/test/univariate/continuous/logitnormal.jl
index 20ccce3bd1..454376606c 100644
--- a/test/univariate/continuous/logitnormal.jl
+++ b/test/univariate/continuous/logitnormal.jl
@@ -43,7 +43,7 @@ function test_logitnormal(g::LogitNormal, n_tsamples::Int=10^6,
     for i = 1:min(100, n_tsamples)
         @test logpdf(g, X[i]) ≈ log(pdf(g, X[i]))
     end
-    @test logpdf.(g, X) ≈ log.(pdf.(g, X))
+    @test Base.Fix1(logpdf, g).(X) ≈ log.(Base.Fix1(pdf, g).(X))
     @test isequal(logpdf(g, 0),-Inf)
     @test isequal(logpdf(g, 1),-Inf)
     @test isequal(logpdf(g, -eps()),-Inf)
diff --git a/test/univariate/continuous/rician.jl b/test/univariate/continuous/rician.jl
index 5674373b36..a75397f891 100644
--- a/test/univariate/continuous/rician.jl
+++ b/test/univariate/continuous/rician.jl
@@ -14,14 +14,14 @@
     @test var(d1) ≈ var(d2)
     @test mode(d1) ≈ mode(d2)
     @test median(d1) ≈ median(d2)
-    @test quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90]) ≈ quantile.(d2, [0.25, 0.45, 0.60, 0.80, 0.90])
-    @test pdf.(d1, 0.0:0.1:1.0) ≈ pdf.(d2, 0.0:0.1:1.0)
-    @test cdf.(d1, 0.0:0.1:1.0) ≈ cdf.(d2, 0.0:0.1:1.0)
+    @test Base.Fix1(quantile, d1).([0.25, 0.45, 0.60, 0.80, 0.90]) ≈ Base.Fix1(quantile, d2).([0.25, 0.45, 0.60, 0.80, 0.90])
+    @test Base.Fix1(pdf, d1).(0.0:0.1:1.0) ≈ Base.Fix1(pdf, d2).(0.0:0.1:1.0)
+    @test Base.Fix1(cdf, d1).(0.0:0.1:1.0) ≈ Base.Fix1(cdf, d2).(0.0:0.1:1.0)
 
     d1 = Rician(10.0, 10.0)
     @test median(d1) == quantile(d1, 0.5)
-    x = quantile.(d1, [0.25, 0.45, 0.60, 0.80, 0.90])
-    @test all(cdf.(d1, x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90])
+    x = Base.Fix1(quantile, d1).([0.25, 0.45, 0.60, 0.80, 0.90])
+    @test all(Base.Fix1(cdf, d1).(x) .≈ [0.25, 0.45, 0.60, 0.80, 0.90])
 
     x = rand(Rician(5.0, 5.0), 100000)
     d1 = fit(Rician, x)
diff --git a/test/univariate/continuous/skewnormal.jl b/test/univariate/continuous/skewnormal.jl
index 81d59ed961..5d1257162c 100644
--- a/test/univariate/continuous/skewnormal.jl
+++ b/test/univariate/continuous/skewnormal.jl
@@ -25,7 +25,7 @@ import Distributions: normpdf, normcdf, normlogpdf, normlogcdf
     d4 = Normal(0.5, 2.2)
     #
     @test pdf(d3, 3.3) == Distributions.pdf(d4, 3.3)
-    @test pdf.(d3, 1:3) == Distributions.pdf.(d4, 1:3)
+    @test Base.Fix1(pdf, d3).(1:3) == Base.Fix1(pdf, d4).(1:3)
     a = mean(d3), var(d3), std(d3)
     b = Distributions.mean(d4), Distributions.var(d4), Distributions.std(d4)
     @test a == b
diff --git a/test/univariate/discrete/binomial.jl b/test/univariate/discrete/binomial.jl
index 9fd3251f2b..0caa974216 100644
--- a/test/univariate/discrete/binomial.jl
+++ b/test/univariate/discrete/binomial.jl
@@ -9,14 +9,14 @@ Random.seed!(1234)
 for (p, n) in [(0.6, 10), (0.8, 6), (0.5, 40), (0.04, 20), (1., 100), (0., 10), (0.999999, 1000), (1e-7, 1000)]
     d = Binomial(n, p)
 
-    a = pdf.(d, 0:n)
+    a = Base.Fix1(pdf, d).(0:n)
     for t=0:n
         @test pdf(d, t) ≈ a[1+t]
     end
 
     li = rand(0:n, 2)
     rng = minimum(li):maximum(li)
-    b = pdf.(d, rng)
+    b = Base.Fix1(pdf, d).(rng)
     for t in rng
         @test pdf(d, t) ≈ b[t - first(rng) + 1]
     end
diff --git a/test/univariate/discrete/categorical.jl b/test/univariate/discrete/categorical.jl
index 45d2c84f7d..a835f7c13c 100644
--- a/test/univariate/discrete/categorical.jl
+++ b/test/univariate/discrete/categorical.jl
@@ -56,8 +56,8 @@ for p in Any[
     @test iszero(ccdf(d, Inf))
     @test isnan(ccdf(d, NaN))
 
-    @test pdf.(d, support(d)) == p
-    @test pdf.(d, 1:k) == p
+    @test Base.Fix1(pdf, d).(support(d)) == p
+    @test Base.Fix1(pdf, d).(1:k) == p
 
     @test cf(d, 0) ≈ 1.0
     @test cf(d, 1) ≈ p' * cis.(1:length(p))
diff --git a/test/univariate/discrete/poissonbinomial.jl b/test/univariate/discrete/poissonbinomial.jl
index 8e53e6edd6..611a54c258 100644
--- a/test/univariate/discrete/poissonbinomial.jl
+++ b/test/univariate/discrete/poissonbinomial.jl
@@ -104,9 +104,9 @@ for (n₁, n₂, n₃, p₁, p₂, p₃) in [(10, 10, 10, 0.1, 0.5, 0.9),
     b2 = Binomial(n₂, p₂)
     b3 = Binomial(n₃, p₃)
 
-    pmf1 = pdf.(b1, support(b1))
-    pmf2 = pdf.(b2, support(b2))
-    pmf3 = pdf.(b3, support(b3))
+    pmf1 = Base.Fix1(pdf, b1).(support(b1))
+    pmf2 = Base.Fix1(pdf, b2).(support(b2))
+    pmf3 = Base.Fix1(pdf, b3).(support(b3))
 
     @test @inferred(mean(d)) ≈ (mean(b1) + mean(b2) + mean(b3))
     @test @inferred(var(d))  ≈ (var(b1) + var(b2) + var(b3))
diff --git a/test/univariate/discrete/soliton.jl b/test/univariate/discrete/soliton.jl
index e0e15e101e..af0aefddc1 100644
--- a/test/univariate/discrete/soliton.jl
+++ b/test/univariate/discrete/soliton.jl
@@ -5,7 +5,7 @@ using Distributions
     Ω = Soliton(K, M, δ, atol)
     @test pdf(Ω, M) > pdf(Ω, M-1)
     @test pdf(Ω, M) > pdf(Ω, M+1)
-    @test cumsum(pdf.(Ω, 1:K)) ≈ cdf.(Ω, 1:K)
+    @test cumsum(Base.Fix1(pdf, Ω).(1:K)) ≈ Base.Fix1(cdf, Ω).(1:K)
     @test cdf(Ω, 0) ≈ 0
     @test cdf(Ω, K) ≈ 1
     @test mean(Ω) ≈ 7.448826535558562
@@ -20,7 +20,7 @@ using Distributions
     K, M, δ, atol = 100, 60, 0.2, 1e-3
     Ω = Soliton(K, M, δ, atol)
     ds = [d for d in 1:K if pdf(Ω, d) > 0]
-    @test all(pdf.(Ω, ds) .> atol)
+    @test all(Base.Fix1(pdf, Ω).(ds) .> atol)
     @test cdf(Ω, 0) ≈ 0
     @test cdf(Ω, K) ≈ 1
     @test minimum(Ω) == 1
diff --git a/test/univariate/locationscale.jl b/test/univariate/locationscale.jl
index 761ce55bbf..c4876e1970 100644
--- a/test/univariate/locationscale.jl
+++ b/test/univariate/locationscale.jl
@@ -83,9 +83,9 @@ function test_location_scale(
             insupport(dtest, -x) == insupport(dref, -x)
 
             @test pdf(dtest, x) ≈ pdf(dref, x)
-            @test pdf.(dtest, xs) ≈ pdf.(dref, xs)
+            @test Base.Fix1(pdf, dtest).(xs) ≈ Base.Fix1(pdf, dref).(xs)
             @test logpdf(dtest, x) ≈ logpdf(dref, x)
-            @test logpdf.(dtest, xs) ≈ logpdf.(dref, xs)
+            @test Base.Fix1(logpdf, dtest).(xs) ≈ Base.Fix1(logpdf, dref).(xs)
             @test loglikelihood(dtest, x) ≈ loglikelihood(dref, x)
             @test loglikelihood(dtest, xs) ≈ loglikelihood(dref, xs)
 
diff --git a/test/univariates.jl b/test/univariates.jl
index e16c52f59b..669ddfe0e7 100644
--- a/test/univariates.jl
+++ b/test/univariates.jl
@@ -99,8 +99,8 @@ function verify_and_test(D::Union{Type,Function}, d::UnivariateDistribution, dct
 
         # pdf method is not implemented for StudentizedRange
         if !isa(d, StudentizedRange)
-            @test isapprox(pdf.(d, x),     p; atol=1e-16, rtol=1e-8)
-            @test isapprox(logpdf.(d, x), lp; atol=isa(d, NoncentralHypergeometric) ? 1e-4 : 1e-12)
+            @test Base.Fix1(pdf, d).(x) ≈ p atol=1e-16 rtol=1e-8
+            @test Base.Fix1(logpdf, d).(x) ≈ lp atol=isa(d, NoncentralHypergeometric) ? 1e-4 : 1e-12
         end
 
         # cdf method is not implemented for NormalInverseGaussian

From b48ca27c9133dfa17d7c384b4e03050e108d041c Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Fri, 5 Jan 2024 11:36:57 +0100
Subject: [PATCH 2/2] Delete src/multivariate/mvnormal copy.jl

---
 src/multivariate/mvnormal copy.jl | 549 ------------------------------
 1 file changed, 549 deletions(-)
 delete mode 100644 src/multivariate/mvnormal copy.jl

diff --git a/src/multivariate/mvnormal copy.jl b/src/multivariate/mvnormal copy.jl
deleted file mode 100644
index 3f0db60f1d..0000000000
--- a/src/multivariate/mvnormal copy.jl	
+++ /dev/null
@@ -1,549 +0,0 @@
-# Multivariate Normal distribution
-
-
-###########################################################
-#
-#   Abstract base class for multivariate normal
-#
-#   Each subtype should provide the following methods:
-#
-#   - length(d):        vector dimension
-#   - mean(d):          the mean vector (in full form)
-#   - cov(d):           the covariance matrix (in full form)
-#   - invcov(d):        inverse of covariance
-#   - logdetcov(d):     log-determinant of covariance
-#   - sqmahal(d, x):        Squared Mahalanobis distance to center
-#   - sqmahal!(r, d, x):    Squared Mahalanobis distances
-#   - gradlogpdf(d, x):     Gradient of logpdf w.r.t. x
-#   - _rand!(d, x):         Sample random vector(s)
-#
-#   Other generic functions will be implemented on top
-#   of these core functions.
-#
-###########################################################
-
-"""
-
-The [Multivariate normal distribution](http://en.wikipedia.org/wiki/Multivariate_normal_distribution)
-is a multidimensional generalization of the *normal distribution*. The probability density function of
-a d-dimensional multivariate normal distribution with mean vector ``\\boldsymbol{\\mu}`` and
-covariance matrix ``\\boldsymbol{\\Sigma}`` is:
-
-```math
-f(\\mathbf{x}; \\boldsymbol{\\mu}, \\boldsymbol{\\Sigma}) = \\frac{1}{(2 \\pi)^{d/2} |\\boldsymbol{\\Sigma}|^{1/2}}
-\\exp \\left( - \\frac{1}{2} (\\mathbf{x} - \\boldsymbol{\\mu})^T \\Sigma^{-1} (\\mathbf{x} - \\boldsymbol{\\mu}) \\right)
-```
-
-We realize that the mean vector and the covariance often have special forms in practice,
-which can be exploited to simplify the computation. For example, the mean vector is sometimes
-just a zero vector, while the covariance matrix can be a diagonal matrix or even in the form
-of ``\\sigma^2 \\mathbf{I}``. To take advantage of such special cases, we introduce a parametric
-type `MvNormal`, defined as below, which allows users to specify the special structure of
-the mean and covariance.
-
-```julia
-struct MvNormal{T<:Real,Cov<:AbstractPDMat,Mean<:AbstractVector} <: AbstractMvNormal
-    μ::Mean
-    Σ::Cov
-end
-```
-
-Here, the mean vector can be an instance of any `AbstractVector`. The covariance can be
-of any subtype of `AbstractPDMat`. Particularly, one can use `PDMat` for full covariance,
-`PDiagMat` for diagonal covariance, and `ScalMat` for the isotropic covariance -- those
-in the form of ``\\sigma^2 \\mathbf{I}``. (See the Julia package
-[PDMats](https://github.com/JuliaStats/PDMats.jl/) for details).
-
-We also define a set of aliases for the types using different combinations of mean vectors and covariance:
-
-```julia
-const IsoNormal  = MvNormal{Float64, ScalMat{Float64},                  Vector{Float64}}
-const DiagNormal = MvNormal{Float64, PDiagMat{Float64,Vector{Float64}}, Vector{Float64}}
-const FullNormal = MvNormal{Float64, PDMat{Float64,Matrix{Float64}},    Vector{Float64}}
-
-const ZeroMeanIsoNormal{Axes}  = MvNormal{Float64, ScalMat{Float64},                  Zeros{Float64,1,Axes}}
-const ZeroMeanDiagNormal{Axes} = MvNormal{Float64, PDiagMat{Float64,Vector{Float64}}, Zeros{Float64,1,Axes}}
-const ZeroMeanFullNormal{Axes} = MvNormal{Float64, PDMat{Float64,Matrix{Float64}},    Zeros{Float64,1,Axes}}
-```
-
-Multivariate normal distributions support affine transformations:
-```julia
-d = MvNormal(μ, Σ)
-c + B * d    # == MvNormal(B * μ + c, B * Σ * B')
-dot(b, d)    # == Normal(dot(b, μ), b' * Σ * b)
-```
-"""
-abstract type AbstractMvNormal <: ContinuousMultivariateDistribution end
-
-### Generic methods (for all AbstractMvNormal subtypes)
-
-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))
-mode(d::AbstractMvNormal) = mean(d)
-modes(d::AbstractMvNormal) = [mean(d)]
-
-"""
-    rand(::AbstractRNG, ::Distributions.AbstractMvNormal)
-
-Sample a random vector from the provided multi-variate normal distribution.
-"""
-rand(::AbstractRNG, ::Distributions.AbstractMvNormal)
-
-function entropy(d::AbstractMvNormal)
-    ldcd = logdetcov(d)
-    return (length(d) * (oftype(ldcd, log2π) + 1) + ldcd) / 2
-end
-
-function mvnormal_c0(d::AbstractMvNormal)
-    ldcd = logdetcov(d)
-    return - (length(d) * oftype(ldcd, log2π) + ldcd) / 2
-end
-
-function kldivergence(p::AbstractMvNormal, q::AbstractMvNormal)
-    # This is the generic implementation for AbstractMvNormal, you might need to specialize for your type
-    length(p) == length(q) ||
-        throw(DimensionMismatch("Distributions p and q have different dimensions $(length(p)) and $(length(q))"))
-    # logdetcov is used separately from _cov for any potential optimization done there
-    return (tr(_cov(q) \ _cov(p)) + sqmahal(q, mean(p)) - length(p) + logdetcov(q) - logdetcov(p)) / 2
-end
-
-# This is a workaround to take advantage of the PDMats objects for MvNormal and avoid copies as Matrix
-# TODO: Remove this once `cov(::MvNormal)` returns the PDMats object
-_cov(d::AbstractMvNormal) = cov(d)
-
-"""
-    invcov(d::AbstractMvNormal)
-
-Return the inversed covariance matrix of d.
-"""
-invcov(d::AbstractMvNormal)
-
-"""
-    logdetcov(d::AbstractMvNormal)
-
-Return the log-determinant value of the covariance matrix.
-"""
-logdetcov(d::AbstractMvNormal)
-
-"""
-    sqmahal(d::AbstractMvNormal, x::AbstractVecOrMat)
-
-Return the squared Mahalanobis distance from `x` to the mean of `d` w.r.t. the covariance.
-When x is a vector, it returns a scalar value. When x is a matrix, it returns a vector of length size(x,2).
-"""
-sqmahal(d::AbstractMvNormal, x::AbstractVecOrMat)
-
-"""
-    sqmahal!(r::AbstractVector, d::AbstractMvNormal, x::AbstractMatrix)
-
-Compute [`sqmahal(d, x)`](@ref) and write the results to the pre-allocated vector `r`.
-"""
-sqmahal!(r::AbstractVector, d::AbstractMvNormal, x::AbstractMatrix)
-
-@inline function logpdf(d::AbstractMvNormal, x::AbstractArray{<:Real,N}) where {N}
-    @boundscheck begin
-        size(x, 1) == length(d) || throw(DimensionMismatch("inconsistent array dimensions"))
-    end
-
-    # Compute constants (with respect to `x`) only once
-    c0 = mvnormal_c0(d) 
-
-    if N == 1
-        # Just numbers in this case
-        return c0 - sqmahal(d, x) / 2
-    else
-        # Compute squared Mahalanobis distance
-        xmat = N == 2 ? x : reshape(x, size(x, 1), :)
-        sqm = sqmahal(d, xmat)
-        
-        # Optimization with reduced allocations for Array
-        res = if sqm isa Array
-            @inbounds for i in eachindex(sqm)
-                sqm[i] = c0 - sqm[i] / 2
-            end
-            sqm
-        else
-            c0 .- sqm ./ 2
-        end
-
-        return N == 2 ? res : reshape(res, ntuple(i -> size(x, i + 1), Val(N - 1)))
-    end
-end
-
-function _logpdf!(r::AbstractArray{<:Real}, d::AbstractMvNormal, x::AbstractArray{<:Real,N}) where {N}
-    dim = length(d)
-    xmat = N == 2 ? x : reshape(x, dim, :)
-    sqmahal!(vec(r), d, xmat)
-    c0 = mvnormal_c0(d)
-    for i in eachindex(r)
-        @inbounds r[i] = c0 - r[i]/2
-    end
-    return r
-end
-
-###########################################################
-#
-#   MvNormal
-#
-#   Multivariate normal distribution with mean parameters
-#
-###########################################################
-"""
-    MvNormal
-
-Generally, users don't have to worry about these internal details.
-
-We provide a common constructor `MvNormal`, which will construct a distribution of
-appropriate type depending on the input arguments.
-"""
-struct MvNormal{T<:Real,Cov<:AbstractPDMat,Mean<:AbstractVector} <: AbstractMvNormal
-    μ::Mean
-    Σ::Cov
-end
-
-const MultivariateNormal = MvNormal  # for the purpose of backward compatibility
-
-const IsoNormal  = MvNormal{Float64,ScalMat{Float64},Vector{Float64}}
-const DiagNormal = MvNormal{Float64,PDiagMat{Float64,Vector{Float64}},Vector{Float64}}
-const FullNormal = MvNormal{Float64,PDMat{Float64,Matrix{Float64}},Vector{Float64}}
-
-const ZeroMeanIsoNormal{Axes}  = MvNormal{Float64,ScalMat{Float64},Zeros{Float64,1,Axes}}
-const ZeroMeanDiagNormal{Axes} = MvNormal{Float64,PDiagMat{Float64,Vector{Float64}},Zeros{Float64,1,Axes}}
-const ZeroMeanFullNormal{Axes} = MvNormal{Float64,PDMat{Float64,Matrix{Float64}},Zeros{Float64,1,Axes}}
-
-### Construction
-function MvNormal(μ::AbstractVector{T}, Σ::AbstractPDMat{T}) where {T<:Real}
-    size(Σ, 1) == length(μ) || throw(DimensionMismatch("The dimensions of mu and Sigma are inconsistent."))
-    MvNormal{T,typeof(Σ), typeof(μ)}(μ, Σ)
-end
-
-function MvNormal(μ::AbstractVector{<:Real}, Σ::AbstractPDMat{<:Real})
-    R = Base.promote_eltype(μ, Σ)
-    MvNormal(convert(AbstractArray{R}, μ), convert(AbstractArray{R}, Σ))
-end
-
-# constructor with general covariance matrix
-"""
-    MvNormal(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real})
-
-Construct a multivariate normal distribution with mean `μ` and covariance matrix `Σ`.
-"""
-MvNormal(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real}) = MvNormal(μ, PDMat(Σ))
-MvNormal(μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real}) = MvNormal(μ, PDiagMat(Σ.diag))
-MvNormal(μ::AbstractVector{<:Real}, Σ::Union{Symmetric{<:Real,<:Diagonal{<:Real}},Hermitian{<:Real,<:Diagonal{<:Real}}}) = MvNormal(μ, PDiagMat(Σ.data.diag))
-MvNormal(μ::AbstractVector{<:Real}, Σ::UniformScaling{<:Real}) =
-    MvNormal(μ, ScalMat(length(μ), Σ.λ))
-function MvNormal(
-    μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real,<:FillArrays.AbstractFill{<:Real,1}}
-)
-    return MvNormal(μ, ScalMat(size(Σ, 1), FillArrays.getindex_value(Σ.diag)))
-end
-
-# constructor without mean vector
-"""
-    MvNormal(Σ::AbstractMatrix{<:Real})
-
-Construct a multivariate normal distribution with zero mean and covariance matrix `Σ`.
-"""
-MvNormal(Σ::AbstractMatrix{<:Real}) = MvNormal(Zeros{eltype(Σ)}(size(Σ, 1)), Σ)
-
-# deprecated constructors with standard deviations
-Base.@deprecate MvNormal(μ::AbstractVector{<:Real}, σ::AbstractVector{<:Real}) MvNormal(μ, LinearAlgebra.Diagonal(map(abs2, σ)))
-Base.@deprecate MvNormal(μ::AbstractVector{<:Real}, σ::Real) MvNormal(μ, σ^2 * I)
-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.Σ))
-end
-Base.convert(::Type{MvNormal{T}}, d::MvNormal{T}) where {T<:Real} = d
-
-function convert(::Type{MvNormal{T}}, μ::AbstractVector, Σ::AbstractPDMat) where T<:Real
-    MvNormal(convert(AbstractArray{T}, μ), convert(AbstractArray{T}, Σ))
-end
-
-### Show
-
-distrname(d::IsoNormal)  = "IsoNormal"    # Note: IsoNormal, etc are just alias names.
-distrname(d::DiagNormal) = "DiagNormal"
-distrname(d::FullNormal) = "FullNormal"
-
-distrname(d::ZeroMeanIsoNormal) = "ZeroMeanIsoNormal"
-distrname(d::ZeroMeanDiagNormal) = "ZeroMeanDiagNormal"
-distrname(d::ZeroMeanFullNormal) = "ZeroMeanFullNormal"
-
-Base.show(io::IO, d::MvNormal) =
-    show_multline(io, d, [(:dim, length(d)), (:μ, mean(d)), (:Σ, cov(d))])
-
-### Basic statistics
-
-length(d::MvNormal) = length(d.μ)
-mean(d::MvNormal) = d.μ
-params(d::MvNormal) = (d.μ, d.Σ)
-@inline partype(d::MvNormal{T}) where {T<:Real} = T
-
-var(d::MvNormal) = diag(d.Σ)
-cov(d::MvNormal) = Matrix(d.Σ)
-_cov(d::MvNormal) = d.Σ
-
-invcov(d::MvNormal) = Matrix(inv(d.Σ))
-logdetcov(d::MvNormal) = logdet(d.Σ)
-
-### Evaluation
-
-sqmahal(d::MvNormal, x::AbstractVector) = invquad(d.Σ, x - d.μ)
-sqmahal(d::MvNormal, x::AbstractMatrix) = invquad(d.Σ, x .- d.μ)
-sqmahal!(r::AbstractVector, d::MvNormal, x::AbstractMatrix) = invquad!(r, d.Σ, x .- d.μ)
-
-gradlogpdf(d::MvNormal, x::AbstractVector{<:Real}) = -(d.Σ \ (x .- d.μ))
-
-# Sampling (for GenericMvNormal)
-
-function _rand!(rng::AbstractRNG, d::MvNormal, x::VecOrMat)
-    unwhiten!(d.Σ, randn!(rng, x))
-    x .+= d.μ
-    return x
-end
-
-# Workaround: randn! only works for Array, but not generally for AbstractArray
-function _rand!(rng::AbstractRNG, d::MvNormal, x::AbstractVector)
-    for i in eachindex(x)
-        @inbounds x[i] = randn(rng, eltype(x))
-    end
-    unwhiten!(d.Σ, x)
-    x .+= d.μ
-    return x
-end
-
-### Affine transformations
-
-Base.:+(d::MvNormal, c::AbstractVector) = MvNormal(d.μ + c, d.Σ)
-Base.:+(c::AbstractVector, d::MvNormal) = d + c
-Base.:-(d::MvNormal, c::AbstractVector) = MvNormal(d.μ - c, d.Σ)
-
-Base.:*(B::AbstractMatrix, d::MvNormal) = MvNormal(B * d.μ, X_A_Xt(d.Σ, B))
-
-dot(b::AbstractVector, d::MvNormal) = Normal(dot(d.μ, b), √quad(d.Σ, b))
-
-dot(d::MvNormal, b::AbstractVector) = dot(b, d)
-
-###########################################################
-#
-#   Estimation of MvNormal
-#
-###########################################################
-
-### Estimation with known covariance
-
-struct MvNormalKnownCov{Cov<:AbstractPDMat}
-    Σ::Cov
-end
-
-MvNormalKnownCov(d::Int, σ::Real) = MvNormalKnownCov(ScalMat(d, abs2(Float64(σ))))
-MvNormalKnownCov(σ::Vector{Float64}) = MvNormalKnownCov(PDiagMat(abs2.(σ)))
-MvNormalKnownCov(Σ::Matrix{Float64}) = MvNormalKnownCov(PDMat(Σ))
-
-length(g::MvNormalKnownCov) = size(g.Σ, 1)
-
-struct MvNormalKnownCovStats{Cov<:AbstractPDMat}
-    invΣ::Cov              # inverse covariance
-    sx::Vector{Float64}    # (weighted) sum of vectors
-    tw::Float64            # sum of weights
-end
-
-function suffstats(g::MvNormalKnownCov{Cov}, x::AbstractMatrix{Float64}) where Cov<:AbstractPDMat
-    size(x,1) == length(g) || throw(DimensionMismatch("Invalid argument dimensions."))
-    invΣ = inv(g.Σ)
-    sx = vec(sum(x, dims=2))
-    tw = Float64(size(x, 2))
-    MvNormalKnownCovStats{Cov}(invΣ, sx, tw)
-end
-
-function suffstats(g::MvNormalKnownCov{Cov}, x::AbstractMatrix{Float64}, w::AbstractVector) where Cov<:AbstractPDMat
-    (size(x,1) == length(g) && size(x,2) == length(w)) ||
-        throw(DimensionMismatch("Inconsistent argument dimensions."))
-    invΣ = inv(g.Σ)
-    sx = x * vec(w)
-    tw = sum(w)
-    MvNormalKnownCovStats{Cov}(invΣ, sx, tw)
-end
-
-## MLE estimation with covariance known
-
-fit_mle(g::MvNormalKnownCov{C}, ss::MvNormalKnownCovStats{C}) where {C<:AbstractPDMat} =
-    MvNormal(ss.sx * inv(ss.tw), g.Σ)
-
-function fit_mle(g::MvNormalKnownCov, x::AbstractMatrix{Float64})
-    d = length(g)
-    size(x,1) == d || throw(DimensionMismatch("Invalid argument dimensions."))
-    μ = lmul!(inv(size(x,2)), vec(sum(x,dims=2)))
-    MvNormal(μ, g.Σ)
-end
-
-function fit_mle(g::MvNormalKnownCov, x::AbstractMatrix{Float64}, w::AbstractVector)
-    d = length(g)
-    (size(x,1) == d && size(x,2) == length(w)) ||
-        throw(DimensionMismatch("Inconsistent argument dimensions."))
-    μ = BLAS.gemv('N', inv(sum(w)), x, vec(w))
-    MvNormal(μ, g.Σ)
-end
-
-
-### Estimation (both mean and cov unknown)
-
-struct MvNormalStats <: SufficientStats
-    s::Vector{Float64}  # (weighted) sum of x
-    m::Vector{Float64}  # (weighted) mean of x
-    s2::Matrix{Float64} # (weighted) sum of (x-μ) * (x-μ)'
-    tw::Float64         # total sample weight
-end
-
-function suffstats(D::Type{MvNormal}, x::AbstractMatrix{Float64})
-    d = size(x, 1)
-    n = size(x, 2)
-    s = vec(sum(x, dims=2))
-    m = s * inv(n)
-    z = x .- m
-    s2 = z * z'
-    MvNormalStats(s, m, s2, Float64(n))
-end
-
-function suffstats(D::Type{MvNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
-    d = size(x, 1)
-    n = size(x, 2)
-    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions."))
-
-    tw = sum(w)
-    s = x * vec(w)
-    m = s * inv(tw)
-    z = similar(x)
-    for j = 1:n
-        xj = view(x,:,j)
-        zj = view(z,:,j)
-        swj = sqrt(w[j])
-        for i = 1:d
-            @inbounds zj[i] = swj * (xj[i] - m[i])
-        end
-    end
-    s2 = z * z'
-    MvNormalStats(s, m, s2, tw)
-end
-
-
-# Maximum Likelihood Estimation
-#
-# Specialized algorithms are respectively implemented for
-# each kind of covariance
-#
-
-fit_mle(D::Type{MvNormal}, ss::MvNormalStats) = fit_mle(FullNormal, ss)
-fit_mle(D::Type{MvNormal}, x::AbstractMatrix{Float64}) = fit_mle(FullNormal, x)
-fit_mle(D::Type{MvNormal}, x::AbstractMatrix{Float64}, w::AbstractArray{Float64}) = fit_mle(FullNormal, x, w)
-
-fit_mle(D::Type{FullNormal}, ss::MvNormalStats) = MvNormal(ss.m, ss.s2 * inv(ss.tw))
-
-function fit_mle(D::Type{FullNormal}, x::AbstractMatrix{Float64})
-    n = size(x, 2)
-    mu = vec(mean(x, dims=2))
-    z = x .- mu
-    C = BLAS.syrk('U', 'N', inv(n), z)
-    LinearAlgebra.copytri!(C, 'U')
-    MvNormal(mu, PDMat(C))
-end
-
-function fit_mle(D::Type{FullNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
-    m = size(x, 1)
-    n = size(x, 2)
-    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions"))
-
-    inv_sw = inv(sum(w))
-    mu = BLAS.gemv('N', inv_sw, x, w)
-
-    z = Matrix{Float64}(undef, m, n)
-    for j = 1:n
-        cj = sqrt(w[j])
-        for i = 1:m
-            @inbounds z[i,j] = (x[i,j] - mu[i]) * cj
-        end
-    end
-    C = BLAS.syrk('U', 'N', inv_sw, z)
-    LinearAlgebra.copytri!(C, 'U')
-    MvNormal(mu, PDMat(C))
-end
-
-function fit_mle(D::Type{DiagNormal}, x::AbstractMatrix{Float64})
-    m = size(x, 1)
-    n = size(x, 2)
-
-    mu = vec(mean(x, dims=2))
-    va = zeros(Float64, m)
-    for j = 1:n
-        for i = 1:m
-            @inbounds va[i] += abs2(x[i,j] - mu[i])
-        end
-    end
-    lmul!(inv(n), va)
-    MvNormal(mu, PDiagMat(va))
-end
-
-function fit_mle(D::Type{DiagNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
-    m = size(x, 1)
-    n = size(x, 2)
-    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions"))
-
-    inv_sw = inv(sum(w))
-    mu = BLAS.gemv('N', inv_sw, x, w)
-
-    va = zeros(Float64, m)
-    for j = 1:n
-        @inbounds wj = w[j]
-        for i = 1:m
-            @inbounds va[i] += abs2(x[i,j] - mu[i]) * wj
-        end
-    end
-    lmul!(inv_sw, va)
-    MvNormal(mu, PDiagMat(va))
-end
-
-function fit_mle(D::Type{IsoNormal}, x::AbstractMatrix{Float64})
-    m = size(x, 1)
-    n = size(x, 2)
-
-    mu = vec(mean(x, dims=2))
-    va = 0.
-    for j = 1:n
-        va_j = 0.
-        for i = 1:m
-            @inbounds va_j += abs2(x[i,j] - mu[i])
-        end
-        va += va_j
-    end
-    MvNormal(mu, ScalMat(m, va / (m * n)))
-end
-
-function fit_mle(D::Type{IsoNormal}, x::AbstractMatrix{Float64}, w::AbstractVector)
-    m = size(x, 1)
-    n = size(x, 2)
-    length(w) == n || throw(DimensionMismatch("Inconsistent argument dimensions"))
-
-    sw = sum(w)
-    inv_sw = 1.0 / sw
-    mu = BLAS.gemv('N', inv_sw, x, w)
-
-    va = 0.
-    for j = 1:n
-        @inbounds wj = w[j]
-        va_j = 0.
-        for i = 1:m
-            @inbounds va_j += abs2(x[i,j] - mu[i]) * wj
-        end
-        va += va_j
-    end
-    MvNormal(mu, ScalMat(m, va / (m * sw)))
-end