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