From d709c34a89020c4bde9b652cc3cb8ef4eabe610f Mon Sep 17 00:00:00 2001 From: Domas Linkevicius Date: Thu, 2 Oct 2025 18:47:43 +0900 Subject: [PATCH 1/8] add mvlaplace --- src/Distributions.jl | 1 + src/multivariate/mvlaplace.jl | 73 ++++++++++++++++++++++++++++++++++ src/multivariates.jl | 3 +- test/multivariate/mvlaplace.jl | 37 +++++++++++++++++ test/runtests.jl | 1 + 5 files changed, 114 insertions(+), 1 deletion(-) create mode 100644 src/multivariate/mvlaplace.jl create mode 100644 test/multivariate/mvlaplace.jl diff --git a/src/Distributions.jl b/src/Distributions.jl index fcf6994347..907480098b 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -139,6 +139,7 @@ export MvNormalCanon, MvNormalKnownCov, MvTDist, + SymmetricMvLaplace, NegativeBinomial, NoncentralBeta, NoncentralChisq, diff --git a/src/multivariate/mvlaplace.jl b/src/multivariate/mvlaplace.jl new file mode 100644 index 0000000000..4a6a8ec53f --- /dev/null +++ b/src/multivariate/mvlaplace.jl @@ -0,0 +1,73 @@ +### Implementation of the symmetric multivariate Laplace distribution using the formatlisation from: +### T. Eltoft, Taesu Kim and Te-Won Lee, +### "On the multivariate Laplace distribution," +### in IEEE Signal Processing Letters, vol. 13, no. 5, pp. 300-303, May 2006, doi: 10.1109/LSP.2006.870353 + + +struct SymmetricMvLaplace{T<:Real,Cov<:AbstractPDMat,iCov<:AbstractMatrix,Mean<:AbstractVector} <: ContinuousMultivariateDistribution + μ::Mean + Γ::Cov + iΓ::Cov + λ::T +end + +### generic methods for SymmetricMvLaplace +length(d::SymmetricMvLaplace) = length(d.μ) +params(d::SymmetricMvLaplace) = (d.μ, d.Σ) + +insupport(d::SymmetricMvLaplace, x::AbstractVector) = + length(d) == length(x) && all(isfinite, x) + +minimum(d::SymmetricMvLaplace) = fill(eltype(d)(-Inf), length(d)) +maximum(d::SymmetricMvLaplace) = fill(eltype(d)(Inf), length(d)) +mode(d::SymmetricMvLaplace) = mean(d) +modes(d::SymmetricMvLaplace) = [mean(d)] +mean(d::SymmetricMvLaplace) = d.μ +var(d::SymmetricMvLaplace) = diag(d.λ * d.Γ) +cov(d::SymmetricMvLaplace) = d.λ * d.Γ +invcov(d::SymmetricMvLaplace) = inv(d.λ * d.Γ) + +### Construction when only an overall covariance is specified +function SymmetricMvLaplace(μ::AbstractVector{T}, Σ::AbstractPDMat{T}) where {T<:Real} + size(Σ, 1) == length(μ) || throw(DimensionMismatch("The dimensions of mu and Sigma are inconsistent.")) + λ = det(Σ)^(1/size(Σ,1)) + Γ = 1/λ * Σ + iΓ = inv(Γ) + SymmetricMvLaplace{T,typeof(Γ), typeof(iΓ), typeof(μ)}(μ, Γ, iΓ, λ) +end + +function SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::AbstractPDMat{<:Real}) + R = Base.promote_eltype(μ, Σ) + SymmetricMvLaplace(convert(AbstractArray{R}, μ), convert(AbstractArray{R}, Σ)) +end + +# constructor with general covariance matrix +SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real}) = SymmetricMvLaplace(μ, PDMat(Σ)) +SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real}) = SymmetricMvLaplace(μ, PDiagMat(Σ.diag)) +SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::Union{Symmetric{<:Real,<:Diagonal{<:Real}},Hermitian{<:Real,<:Diagonal{<:Real}}}) = SymmetricMvLaplace(μ, PDiagMat(Σ.data.diag)) +SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::UniformScaling{<:Real}) = + SymmetricMvLaplace(μ, ScalMat(length(μ), Σ.λ)) +function SymmetricMvLaplace( + μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real,<:FillArrays.AbstractFill{<:Real,1}} +) + return SymmetricMvLaplace(μ, ScalMat(size(Σ, 1), FillArrays.getindex_value(Σ.diag))) +end + +### generic methods for SymmetricMvLaplace +Base.eltype(::Type{<:SymmetricMvLaplace{T}}) where {T} = T +@inline partype(d::SymmetricMvLaplace{T}) where {T<:Real} = T + +function _rand!(rng::AbstractRNG, d::SymmetricMvLaplace, x::VecOrMat) + unwhiten!(d.Γ, randn!(rng, x)) + x .*= sqrt.(rand(Exponential(d.λ), 1, size(x,2))) + x .+= d.μ + return x +end + +function _logpdf(d::SymmetricMvLaplace, x::AbstractArray) + dim_half = length(d) / 2 + _2byλ = 2/d.λ + xdif = x - d.μ + _qs = sqrt(_2byλ * dot(xdif, d.iΓ, xdif)) + return log(2π^(-dim_half)) + log(_2byλ) + bessely(dim_half-1, _qs) - log(_qs^(dim_half - 1)) +end \ No newline at end of file diff --git a/src/multivariates.jl b/src/multivariates.jl index 4d98ebe110..7fabb24010 100644 --- a/src/multivariates.jl +++ b/src/multivariates.jl @@ -126,6 +126,7 @@ for fname in ["dirichlet.jl", "mvlognormal.jl", "mvtdist.jl", "product.jl", # deprecated - "vonmisesfisher.jl"] + "vonmisesfisher.jl", + "mvlaplace.jl"] include(joinpath("multivariate", fname)) end diff --git a/test/multivariate/mvlaplace.jl b/test/multivariate/mvlaplace.jl new file mode 100644 index 0000000000..02e3b021b4 --- /dev/null +++ b/test/multivariate/mvlaplace.jl @@ -0,0 +1,37 @@ +using Distributions +import PDMats: ScalMat, PDiagMat, PDMat +using LinearAlgebra, Test +using FillArrays + + +@testset "SymmetricMvLaplace tests" begin + mu = [1., 2., 3.] + mu_r = 1.:3. + + C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.] + dv = [1.2, 3.4, 2.6] + for (g, μ, Σ) in [ + (SymmetricMvLaplace(mu, C), mu, C), + (SymmetricMvLaplace(mu_r, C), mu_r, C), + (SymmetricMvLaplace(mu, Symmetric(C)), mu, Matrix(Symmetric(C))), + (SymmetricMvLaplace(mu_r, Symmetric(C)), mu_r, Matrix(Symmetric(C))), + (SymmetricMvLaplace(mu, Diagonal(dv)), mu, Matrix(Diagonal(dv))), + (SymmetricMvLaplace(mu, Symmetric(Diagonal(dv))), mu, Matrix(Diagonal(dv))), + (SymmetricMvLaplace(mu, Hermitian(Diagonal(dv))), mu, Matrix(Diagonal(dv))), + (SymmetricMvLaplace(mu_r, Diagonal(dv)), mu_r, Matrix(Diagonal(dv))) ] + + @test mean(g) ≈ μ + @test cov(g) ≈ Σ + @test invcov(g) ≈ inv(Σ) + end +end + +@testset "SymmetricMvLaplace constructor" begin + mu = [1., 2., 3.] + C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.] + @test typeof(SymmetricMvLaplace(mu, PDMat(Array{Float32}(C)))) == typeof(SymmetricMvLaplace(mu, PDMat(C))) + @test typeof(SymmetricMvLaplace(mu, Array{Float32}(C))) == typeof(SymmetricMvLaplace(mu, PDMat(C))) + @test SymmetricMvLaplace(mu, I) === SymmetricMvLaplace(mu, Diagonal(Ones(length(mu)))) + @test SymmetricMvLaplace(mu, 9 * I) === SymmetricMvLaplace(mu, Diagonal(Fill(9, length(mu)))) + @test SymmetricMvLaplace(mu, 0.25f0 * I) === SymmetricMvLaplace(mu, Diagonal(Fill(0.25f0, length(mu)))) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index f63a3d5e8e..b52e31c6d1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -46,6 +46,7 @@ const tests = [ "multivariate/dirichletmultinomial", "univariate/continuous/logitnormal", "multivariate/mvtdist", + "multivariate/mvlaplace", "univariate/continuous/kolmogorov", "edgeworth", "matrixreshaped", # extra file compared to /src From cb16f3698f53fcf5eed8deac63e513b8e1c0b0a7 Mon Sep 17 00:00:00 2001 From: Domas Linkevicius Date: Thu, 2 Oct 2025 19:42:36 +0900 Subject: [PATCH 2/8] add docstring --- src/multivariate/mvlaplace.jl | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/multivariate/mvlaplace.jl b/src/multivariate/mvlaplace.jl index 4a6a8ec53f..faeaf55de1 100644 --- a/src/multivariate/mvlaplace.jl +++ b/src/multivariate/mvlaplace.jl @@ -2,8 +2,35 @@ ### T. Eltoft, Taesu Kim and Te-Won Lee, ### "On the multivariate Laplace distribution," ### in IEEE Signal Processing Letters, vol. 13, no. 5, pp. 300-303, May 2006, doi: 10.1109/LSP.2006.870353 +""" +The [symmetric multivariate Laplace distribution](https://en.wikipedia.org/wiki/Multivariate_Laplace_distribution#Symmetric_multivariate_Laplace_distribution) +is a multidimensional generalization of the *Laplace distribution*, as described in +T. Eltoft, Taesu Kim and Te-Won Lee, +"On the multivariate Laplace distribution," +in IEEE Signal Processing Letters, vol. 13, no. 5, pp. 300-303, May 2006, doi: 10.1109/LSP.2006.870353. +The symmetric multivariate Laplace distribution is defined by three parameters: +- ``\lambda``, which is a positive real number used to define an exponential distribution `Exp(λ)` +- ``\\boldsymbol{\\Gamma}``, which is a k-by-k positive definite matrix with `det(Γ)=1` (as per the assumptions in the source paper) +- ``\\boldsymbol{\\mu}``, which is a k-dimensional real-valued vector +The expected valued of the symmetric multivariate Laplace distribution is simply ``\\boldsymbol{\\mu}``, +whereas the covariance ``\\boldsymbol{\\Sigma} = \lambda \\boldsymbol{\\Gamma}``. +The symmetric multivariate Laplace distribution can be created by specifying either a ``\\boldsymbol{\\mu}`` and a ``\\boldsymbol{\\Sigma}`` +(analogously to a `MvNormal`) and the ``\lambda`` and ``\\boldsymbol{\\Gamma}`` are calculated internally, +or by specifying all three parameters, ``\\boldsymbol{\\mu}``, ``\lambda`` and ``\\boldsymbol{\\Gamma}``. + +The probability density function of +a k-dimensional symmetric multivariate Laplace distribution with parameters ``\\boldsymbol{\\mu}``, ``\lambda`` and ``\\boldsymbol{\\Gamma}`` is: +```math +f(\\mathbf{x}; \\boldsymbol{\\mu}, \\lambda, \\boldsymbol{\\Gamma}) = \\frac{1}{(2 \\pi)^{d/2}} \\frac{2}{\\lambda} +\frac{\\mathrm{K}_{(d/2)-1}\\left(\\sqrt{\\frac{2}{\\lambda}q(\\mathbf{x})}\\right)}{\\left(\\sqrt{\\frac{2}{\\lambda}q(\\mathbf{x})}\\right)^{(d/2)-1}} +``` +where ``\\mathrm{K}_m (y)`` is the Bessel function of the second kind of order ``m`` evaluated at ``y`` and +```math +q(\\mathbf{x} = (\\mathbf{x} - \\boldsymbol{\\mu})^T \\Gamma^{-1} (\\mathbf{x} - \\boldsymbol{\\mu}) +``` +""" struct SymmetricMvLaplace{T<:Real,Cov<:AbstractPDMat,iCov<:AbstractMatrix,Mean<:AbstractVector} <: ContinuousMultivariateDistribution μ::Mean Γ::Cov From 324e1d67ab9ebd7c8c275cb3fe33277c85bbf450 Mon Sep 17 00:00:00 2001 From: Domas Linkevicius Date: Thu, 2 Oct 2025 20:03:21 +0900 Subject: [PATCH 3/8] add extra constructor --- src/multivariate/mvlaplace.jl | 38 ++++++++++++++++++++++++++--------- 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/src/multivariate/mvlaplace.jl b/src/multivariate/mvlaplace.jl index faeaf55de1..72c61ad617 100644 --- a/src/multivariate/mvlaplace.jl +++ b/src/multivariate/mvlaplace.jl @@ -1,7 +1,3 @@ -### Implementation of the symmetric multivariate Laplace distribution using the formatlisation from: -### T. Eltoft, Taesu Kim and Te-Won Lee, -### "On the multivariate Laplace distribution," -### in IEEE Signal Processing Letters, vol. 13, no. 5, pp. 300-303, May 2006, doi: 10.1109/LSP.2006.870353 """ The [symmetric multivariate Laplace distribution](https://en.wikipedia.org/wiki/Multivariate_Laplace_distribution#Symmetric_multivariate_Laplace_distribution) @@ -10,18 +6,18 @@ T. Eltoft, Taesu Kim and Te-Won Lee, "On the multivariate Laplace distribution," in IEEE Signal Processing Letters, vol. 13, no. 5, pp. 300-303, May 2006, doi: 10.1109/LSP.2006.870353. The symmetric multivariate Laplace distribution is defined by three parameters: -- ``\lambda``, which is a positive real number used to define an exponential distribution `Exp(λ)` +- ``\\lambda``, which is a positive real number used to define an exponential distribution `Exp(λ)` - ``\\boldsymbol{\\Gamma}``, which is a k-by-k positive definite matrix with `det(Γ)=1` (as per the assumptions in the source paper) - ``\\boldsymbol{\\mu}``, which is a k-dimensional real-valued vector The expected valued of the symmetric multivariate Laplace distribution is simply ``\\boldsymbol{\\mu}``, -whereas the covariance ``\\boldsymbol{\\Sigma} = \lambda \\boldsymbol{\\Gamma}``. +whereas the covariance ``\\boldsymbol{\\Sigma} = \\lambda \\boldsymbol{\\Gamma}``. The symmetric multivariate Laplace distribution can be created by specifying either a ``\\boldsymbol{\\mu}`` and a ``\\boldsymbol{\\Sigma}`` -(analogously to a `MvNormal`) and the ``\lambda`` and ``\\boldsymbol{\\Gamma}`` are calculated internally, -or by specifying all three parameters, ``\\boldsymbol{\\mu}``, ``\lambda`` and ``\\boldsymbol{\\Gamma}``. +(analogously to a `MvNormal`) and the ``\\lambda`` and ``\\boldsymbol{\\Gamma}`` are calculated internally, +or by specifying all three parameters, ``\\boldsymbol{\\mu}``, ``\\lambda`` and ``\\boldsymbol{\\Gamma}``. The probability density function of -a k-dimensional symmetric multivariate Laplace distribution with parameters ``\\boldsymbol{\\mu}``, ``\lambda`` and ``\\boldsymbol{\\Gamma}`` is: +a k-dimensional symmetric multivariate Laplace distribution with parameters ``\\boldsymbol{\\mu}``, ``\\lambda`` and ``\\boldsymbol{\\Gamma}`` is: ```math f(\\mathbf{x}; \\boldsymbol{\\mu}, \\lambda, \\boldsymbol{\\Gamma}) = \\frac{1}{(2 \\pi)^{d/2}} \\frac{2}{\\lambda} \frac{\\mathrm{K}_{(d/2)-1}\\left(\\sqrt{\\frac{2}{\\lambda}q(\\mathbf{x})}\\right)}{\\left(\\sqrt{\\frac{2}{\\lambda}q(\\mathbf{x})}\\right)^{(d/2)-1}} @@ -54,6 +50,18 @@ var(d::SymmetricMvLaplace) = diag(d.λ * d.Γ) cov(d::SymmetricMvLaplace) = d.λ * d.Γ invcov(d::SymmetricMvLaplace) = inv(d.λ * d.Γ) +### Construction when all three parameters are specified +function SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::AbstractPDMat{<:Real}) + size(Γ, 1) == length(μ) || throw(DimensionMismatch("The dimensions of mu and Gamma are inconsistent.")) + isapprox(det(Γ), 1) || throw(ArgumentError("det(Gamma) is not approximately equal to 1.")) + λ > 0 || throw(ArgumentError("λ is not a positive real number.")) + + R = Base.promote_eltype(μ, λ, Γ) + _μ, _λ, _Γ = convert(AbstractArray{R}, μ), convert(R, λ), convert(AbstractArray{R}, Σ) + _iΓ = inv(_Γ) + SymmetricMvLaplace{R,typeof(_Γ), typeof(_iΓ), typeof(_μ)}(_μ, _Γ, _iΓ, _λ) +end + ### Construction when only an overall covariance is specified function SymmetricMvLaplace(μ::AbstractVector{T}, Σ::AbstractPDMat{T}) where {T<:Real} size(Σ, 1) == length(μ) || throw(DimensionMismatch("The dimensions of mu and Sigma are inconsistent.")) @@ -68,6 +76,18 @@ function SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::AbstractPDMat{<:Real SymmetricMvLaplace(convert(AbstractArray{R}, μ), convert(AbstractArray{R}, Σ)) end +# constructor with general gamma matrix and lambda +SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::AbstractMatrix{<:Real}) = SymmetricMvLaplace(μ, λ, PDMat(Γ)) +SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::Diagonal{<:Real}) = SymmetricMvLaplace(μ, λ, PDiagMat(Γ.diag)) +SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::Union{Symmetric{<:Real,<:Diagonal{<:Real}},Hermitian{<:Real,<:Diagonal{<:Real}}}) = SymmetricMvLaplace(μ, λ, PDiagMat(Γ.data.diag)) +SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::UniformScaling{<:Real}) = + SymmetricMvLaplace(μ, λ, ScalMat(length(μ), Γ.λ)) +function SymmetricMvLaplace( + μ::AbstractVector{<:Real}, λ::Real, Γ::Diagonal{<:Real,<:FillArrays.AbstractFill{<:Real,1}} +) + return SymmetricMvLaplace(μ, λ, ScalMat(size(Γ, 1), FillArrays.getindex_value(Γ.diag))) +end + # constructor with general covariance matrix SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real}) = SymmetricMvLaplace(μ, PDMat(Σ)) SymmetricMvLaplace(μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real}) = SymmetricMvLaplace(μ, PDiagMat(Σ.diag)) From 52e8a0868b8c7eeadd7ff7ba4223534f64560451 Mon Sep 17 00:00:00 2001 From: Domas Linkevicius Date: Fri, 3 Oct 2025 17:01:26 +0900 Subject: [PATCH 4/8] expand tests --- src/multivariate/mvlaplace.jl | 13 +++---- test/multivariate/mvlaplace.jl | 68 +++++++++++++++++++++++++++++----- 2 files changed, 65 insertions(+), 16 deletions(-) diff --git a/src/multivariate/mvlaplace.jl b/src/multivariate/mvlaplace.jl index 72c61ad617..009a877d54 100644 --- a/src/multivariate/mvlaplace.jl +++ b/src/multivariate/mvlaplace.jl @@ -36,7 +36,7 @@ end ### generic methods for SymmetricMvLaplace length(d::SymmetricMvLaplace) = length(d.μ) -params(d::SymmetricMvLaplace) = (d.μ, d.Σ) +params(d::SymmetricMvLaplace) = (d.μ, d.λ, d.Γ) insupport(d::SymmetricMvLaplace, x::AbstractVector) = length(d) == length(x) && all(isfinite, x) @@ -57,7 +57,7 @@ function SymmetricMvLaplace(μ::AbstractVector{<:Real}, λ::Real, Γ::AbstractPD λ > 0 || throw(ArgumentError("λ is not a positive real number.")) R = Base.promote_eltype(μ, λ, Γ) - _μ, _λ, _Γ = convert(AbstractArray{R}, μ), convert(R, λ), convert(AbstractArray{R}, Σ) + _μ, _λ, _Γ = convert(AbstractArray{R}, μ), convert(R, λ), convert(AbstractArray{R}, Γ) _iΓ = inv(_Γ) SymmetricMvLaplace{R,typeof(_Γ), typeof(_iΓ), typeof(_μ)}(_μ, _Γ, _iΓ, _λ) end @@ -106,15 +106,14 @@ Base.eltype(::Type{<:SymmetricMvLaplace{T}}) where {T} = T function _rand!(rng::AbstractRNG, d::SymmetricMvLaplace, x::VecOrMat) unwhiten!(d.Γ, randn!(rng, x)) - x .*= sqrt.(rand(Exponential(d.λ), 1, size(x,2))) + x .*= sqrt.(rand(rng, Exponential(d.λ), 1, size(x,2))) x .+= d.μ return x end function _logpdf(d::SymmetricMvLaplace, x::AbstractArray) - dim_half = length(d) / 2 - _2byλ = 2/d.λ + _d = length(d) / 2 xdif = x - d.μ - _qs = sqrt(_2byλ * dot(xdif, d.iΓ, xdif)) - return log(2π^(-dim_half)) + log(_2byλ) + bessely(dim_half-1, _qs) - log(_qs^(dim_half - 1)) + q = dot(xdif, d.iΓ, xdif) + return _d * log(2π) + log(2/d.λ) + log(besselk(_d-1, sqrt(2/d.λ * q))) + 0.5*(_d-1)*log(0.5*d.λ*q) end \ No newline at end of file diff --git a/test/multivariate/mvlaplace.jl b/test/multivariate/mvlaplace.jl index 02e3b021b4..b381824598 100644 --- a/test/multivariate/mvlaplace.jl +++ b/test/multivariate/mvlaplace.jl @@ -26,12 +26,62 @@ using FillArrays end end -@testset "SymmetricMvLaplace constructor" begin - mu = [1., 2., 3.] - C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.] - @test typeof(SymmetricMvLaplace(mu, PDMat(Array{Float32}(C)))) == typeof(SymmetricMvLaplace(mu, PDMat(C))) - @test typeof(SymmetricMvLaplace(mu, Array{Float32}(C))) == typeof(SymmetricMvLaplace(mu, PDMat(C))) - @test SymmetricMvLaplace(mu, I) === SymmetricMvLaplace(mu, Diagonal(Ones(length(mu)))) - @test SymmetricMvLaplace(mu, 9 * I) === SymmetricMvLaplace(mu, Diagonal(Fill(9, length(mu)))) - @test SymmetricMvLaplace(mu, 0.25f0 * I) === SymmetricMvLaplace(mu, Diagonal(Fill(0.25f0, length(mu)))) -end \ No newline at end of file +@testset "SymmetricMvLaplace constructors" begin + @testset "Providing mu and Sigma" begin + mu = [1., 2., 3.] + C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.] + @test typeof(SymmetricMvLaplace(mu, PDMat(Array{Float32}(C)))) == typeof(SymmetricMvLaplace(mu, PDMat(C))) + @test typeof(SymmetricMvLaplace(mu, Array{Float32}(C))) == typeof(SymmetricMvLaplace(mu, PDMat(C))) + @test SymmetricMvLaplace(mu, I) === SymmetricMvLaplace(mu, Diagonal(Ones(length(mu)))) + @test SymmetricMvLaplace(mu, 9 * I) === SymmetricMvLaplace(mu, Diagonal(Fill(9, length(mu)))) + @test SymmetricMvLaplace(mu, 0.25f0 * I) === SymmetricMvLaplace(mu, Diagonal(Fill(0.25f0, length(mu)))) + end + @testset "Providing mu, lambda and Gamma" begin + mu = [1., 2., 3.] + C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.] + L = det(C)^(1/size(C,1)) + G = 1/L * C + @test typeof(SymmetricMvLaplace(mu, L, PDMat(Array{Float32}(G)))) == typeof(SymmetricMvLaplace(mu, L, PDMat(G))) + @test typeof(SymmetricMvLaplace(mu, L, Array{Float32}(G))) == typeof(SymmetricMvLaplace(mu, L, PDMat(G))) + @test SymmetricMvLaplace(mu, 1, I) === SymmetricMvLaplace(mu, Diagonal(Ones(length(mu)))) + + L9 = det(9*I(3))^(1/3) + L025 = det(0.25f0*I(3))^(1/3) + @test SymmetricMvLaplace(mu, L9, I) === SymmetricMvLaplace(mu, L9, Diagonal(Fill(1, length(mu)))) + @test SymmetricMvLaplace(mu, L025, I) === SymmetricMvLaplace(mu, L025, Diagonal(Fill(1, length(mu)))) + end +end + +@testset "Symmetric MvLaplace basic functions" begin + for t in (Float32, Float64) + @testset "Type: $t" begin + mu = t.(ones(3)) + λ = t.(1) + Γ = t.(I(3)) + d = SymmetricMvLaplace(mu, λ, Γ) + + @test length(d) == 3 + @test params(d) == (mu, λ, Γ) + @test maximum(d) == [Inf, Inf, Inf] + @test minimum(d) == [-Inf, -Inf, -Inf] + @test mode(d) == mu + @test modes(d) == [mu] + @test var(d) == diag(λ*Γ) + @test eltype(d) == t + @test partype(d) == t + + s = rand(d, 10) + @test size(s) == (3, 10) + @test eltype(s) == t + @test isfinite(loglikelihood(d, s)) + + l = Laplace(1, 0.707) + s_d = rand(d, 10_000_000) + qs = collect(0.01:0.01:0.99) + q_d = quantile.(Ref(s_d[1,:]), qs) + q_l = quantile.(Ref(l), qs) + @test all(isapprox.(q_d, q_l; atol=1e-2)) ### tolerance a bit high as tails are heavy and a bit stochastic + ### probably need to come up with some logpdf correctness tests too + end + end +end From a1da067c33bb4e30ae9b28d36bb92a0898ded978 Mon Sep 17 00:00:00 2001 From: Domas Linkevicius Date: Sat, 4 Oct 2025 11:05:05 +0900 Subject: [PATCH 5/8] tests to cover the final few lines --- test/multivariate/mvlaplace.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/multivariate/mvlaplace.jl b/test/multivariate/mvlaplace.jl index b381824598..2995561a27 100644 --- a/test/multivariate/mvlaplace.jl +++ b/test/multivariate/mvlaplace.jl @@ -39,6 +39,7 @@ end @testset "Providing mu, lambda and Gamma" begin mu = [1., 2., 3.] C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.] + Csym = C'*C L = det(C)^(1/size(C,1)) G = 1/L * C @test typeof(SymmetricMvLaplace(mu, L, PDMat(Array{Float32}(G)))) == typeof(SymmetricMvLaplace(mu, L, PDMat(G))) @@ -47,8 +48,10 @@ end L9 = det(9*I(3))^(1/3) L025 = det(0.25f0*I(3))^(1/3) + LC = det(Csym)^(1/3) @test SymmetricMvLaplace(mu, L9, I) === SymmetricMvLaplace(mu, L9, Diagonal(Fill(1, length(mu)))) @test SymmetricMvLaplace(mu, L025, I) === SymmetricMvLaplace(mu, L025, Diagonal(Fill(1, length(mu)))) + @test SymmetricMvLaplace(mu, LC, 1/LC*Symmetric(Csym)) == SymmetricMvLaplace(mu, LC, 1/LC * Csym) end end @@ -69,6 +72,7 @@ end @test var(d) == diag(λ*Γ) @test eltype(d) == t @test partype(d) == t + @test insupport(d, rand(3)) s = rand(d, 10) @test size(s) == (3, 10) From 9af30750a7ba2bb32198eacd758dd1c2cc2af47a Mon Sep 17 00:00:00 2001 From: Domas Linkevicius Date: Sat, 4 Oct 2025 12:40:36 +0900 Subject: [PATCH 6/8] adjust test --- test/multivariate/mvlaplace.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/multivariate/mvlaplace.jl b/test/multivariate/mvlaplace.jl index 2995561a27..d54c0fd74d 100644 --- a/test/multivariate/mvlaplace.jl +++ b/test/multivariate/mvlaplace.jl @@ -39,7 +39,6 @@ end @testset "Providing mu, lambda and Gamma" begin mu = [1., 2., 3.] C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.] - Csym = C'*C L = det(C)^(1/size(C,1)) G = 1/L * C @test typeof(SymmetricMvLaplace(mu, L, PDMat(Array{Float32}(G)))) == typeof(SymmetricMvLaplace(mu, L, PDMat(G))) @@ -48,10 +47,9 @@ end L9 = det(9*I(3))^(1/3) L025 = det(0.25f0*I(3))^(1/3) - LC = det(Csym)^(1/3) @test SymmetricMvLaplace(mu, L9, I) === SymmetricMvLaplace(mu, L9, Diagonal(Fill(1, length(mu)))) @test SymmetricMvLaplace(mu, L025, I) === SymmetricMvLaplace(mu, L025, Diagonal(Fill(1, length(mu)))) - @test SymmetricMvLaplace(mu, LC, 1/LC*Symmetric(Csym)) == SymmetricMvLaplace(mu, LC, 1/LC * Csym) + @test SymmetricMvLaplace(mu, L9, Symmetric(Float64.(I(3)))) == SymmetricMvLaplace(mu, L9, Float64.(I(3))) end end From 09e5bfaf14a997c560bbc7660ce7e9233492e9b4 Mon Sep 17 00:00:00 2001 From: Domas Linkevicius Date: Sun, 5 Oct 2025 07:23:24 +0900 Subject: [PATCH 7/8] add logpdf correctness tests --- src/multivariate/mvlaplace.jl | 5 ++++- test/multivariate/mvlaplace.jl | 23 +++++++++++++---------- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/src/multivariate/mvlaplace.jl b/src/multivariate/mvlaplace.jl index 009a877d54..89f099ca42 100644 --- a/src/multivariate/mvlaplace.jl +++ b/src/multivariate/mvlaplace.jl @@ -26,6 +26,9 @@ where ``\\mathrm{K}_m (y)`` is the Bessel function of the second kind of order ` ```math q(\\mathbf{x} = (\\mathbf{x} - \\boldsymbol{\\mu})^T \\Gamma^{-1} (\\mathbf{x} - \\boldsymbol{\\mu}) ``` +Note that the density function becomes undefined as ``x\\rightarrow\mu`` because ``\\mathrm{K}_m (y)`` approaches +infinity as ``y`` approaches zero. This feature is not unique to the formulation from Eltoft et al., but is also present in +the (slightly different) version presented in https://en.wikipedia.org/wiki/Multivariate_Laplace_distribution. """ struct SymmetricMvLaplace{T<:Real,Cov<:AbstractPDMat,iCov<:AbstractMatrix,Mean<:AbstractVector} <: ContinuousMultivariateDistribution μ::Mean @@ -115,5 +118,5 @@ function _logpdf(d::SymmetricMvLaplace, x::AbstractArray) _d = length(d) / 2 xdif = x - d.μ q = dot(xdif, d.iΓ, xdif) - return _d * log(2π) + log(2/d.λ) + log(besselk(_d-1, sqrt(2/d.λ * q))) + 0.5*(_d-1)*log(0.5*d.λ*q) + return -_d * log(2π) + log(2/d.λ) + log(besselk(_d-1, sqrt(2/d.λ * q))) - 0.5*(_d-1)*log(0.5*d.λ*q) end \ No newline at end of file diff --git a/test/multivariate/mvlaplace.jl b/test/multivariate/mvlaplace.jl index d54c0fd74d..5b928bbe97 100644 --- a/test/multivariate/mvlaplace.jl +++ b/test/multivariate/mvlaplace.jl @@ -54,36 +54,39 @@ end end @testset "Symmetric MvLaplace basic functions" begin + d_dim = 2 for t in (Float32, Float64) @testset "Type: $t" begin - mu = t.(ones(3)) + mu = t.(zeros(d_dim)) λ = t.(1) - Γ = t.(I(3)) + Γ = t.(I(d_dim)) d = SymmetricMvLaplace(mu, λ, Γ) - @test length(d) == 3 + @test length(d) == d_dim @test params(d) == (mu, λ, Γ) - @test maximum(d) == [Inf, Inf, Inf] - @test minimum(d) == [-Inf, -Inf, -Inf] + @test maximum(d) == repeat([Inf], d_dim) + @test minimum(d) == repeat([-Inf], d_dim) @test mode(d) == mu @test modes(d) == [mu] @test var(d) == diag(λ*Γ) @test eltype(d) == t @test partype(d) == t - @test insupport(d, rand(3)) + @test insupport(d, rand(2)) s = rand(d, 10) - @test size(s) == (3, 10) + @test size(s) == (2, 10) @test eltype(s) == t @test isfinite(loglikelihood(d, s)) - l = Laplace(1, 0.707) + l = Laplace(0, 0.707) s_d = rand(d, 10_000_000) qs = collect(0.01:0.01:0.99) q_d = quantile.(Ref(s_d[1,:]), qs) q_l = quantile.(Ref(l), qs) @test all(isapprox.(q_d, q_l; atol=1e-2)) ### tolerance a bit high as tails are heavy and a bit stochastic - ### probably need to come up with some logpdf correctness tests too + ### below rhs calculated by pen and paper + @test isapprox(logpdf(d, [1; 1]), log(1/π * besselk(0, 2))) + @test isapprox(logpdf(d, [2; 2]), log(1/π * besselk(0, 4))) end end -end +end \ No newline at end of file From b3a6400acb09463a78868ef013b46a751487a47b Mon Sep 17 00:00:00 2001 From: Domas Linkevicius Date: Sun, 5 Oct 2025 07:26:25 +0900 Subject: [PATCH 8/8] fix typo --- src/multivariate/mvlaplace.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/multivariate/mvlaplace.jl b/src/multivariate/mvlaplace.jl index 89f099ca42..9ba8a17c2e 100644 --- a/src/multivariate/mvlaplace.jl +++ b/src/multivariate/mvlaplace.jl @@ -26,7 +26,7 @@ where ``\\mathrm{K}_m (y)`` is the Bessel function of the second kind of order ` ```math q(\\mathbf{x} = (\\mathbf{x} - \\boldsymbol{\\mu})^T \\Gamma^{-1} (\\mathbf{x} - \\boldsymbol{\\mu}) ``` -Note that the density function becomes undefined as ``x\\rightarrow\mu`` because ``\\mathrm{K}_m (y)`` approaches +Note that the density function becomes undefined as ``x\\rightarrow\\boldsymbol{\\mu}`` because ``\\mathrm{K}_m (y)`` approaches infinity as ``y`` approaches zero. This feature is not unique to the formulation from Eltoft et al., but is also present in the (slightly different) version presented in https://en.wikipedia.org/wiki/Multivariate_Laplace_distribution. """