From d285de5e0df4290905c0eab5ed0c94ecd1037567 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mos=C3=A8=20Giordano?= Date: Thu, 27 May 2021 11:57:06 +0100 Subject: [PATCH] Add function to compute L2 norm of distributions (#1321) Co-authored-by: David Widmann --- Project.toml | 2 +- docs/src/univariate.md | 1 + src/Distributions.jl | 2 ++ src/pdfnorm.jl | 50 +++++++++++++++++++++++++++++++++ test/pdfnorm.jl | 63 ++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 6 files changed, 118 insertions(+), 1 deletion(-) create mode 100644 src/pdfnorm.jl create mode 100644 test/pdfnorm.jl diff --git a/Project.toml b/Project.toml index 244d2738b..c5b7c37af 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.1" +version = "0.25.2" [deps] FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" diff --git a/docs/src/univariate.md b/docs/src/univariate.md index 61f66dd7d..0418a38fa 100644 --- a/docs/src/univariate.md +++ b/docs/src/univariate.md @@ -63,6 +63,7 @@ entropy(::UnivariateDistribution, ::Bool) entropy(::UnivariateDistribution, ::Real) mgf(::UnivariateDistribution, ::Any) cf(::UnivariateDistribution, ::Any) +pdfsquaredL2norm ``` ### Probability Evaluation diff --git a/src/Distributions.jl b/src/Distributions.jl index 9196e8240..72005e7f8 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -230,6 +230,7 @@ export params!, # provide storage space to calculate the tuple of parameters for a multivariate distribution like mvlognormal partype, # returns a type large enough to hold all of a distribution's parameters' element types pdf, # probability density function (ContinuousDistribution) + pdfsquaredL2norm, # squared L2 norm of the probability density function probs, # Get the vector of probabilities probval, # The pdf/pmf value for a uniform distribution product_distribution, # product of univariate distributions @@ -285,6 +286,7 @@ include("conversion.jl") include("convolution.jl") include("qq.jl") include("estimators.jl") +include("pdfnorm.jl") # mixture distributions (TODO: moveout) include("mixtures/mixturemodel.jl") diff --git a/src/pdfnorm.jl b/src/pdfnorm.jl new file mode 100644 index 000000000..28abcabfd --- /dev/null +++ b/src/pdfnorm.jl @@ -0,0 +1,50 @@ +""" + pdfsquaredL2norm(d::Distribution) + +Return the square of the L2 norm of the probability density function ``f(x)`` of the distribution `d`: + +```math +\\int_{S} f(x)^{2} \\mathrm{d} x +``` + +where ``S`` is the support of ``f(x)``. +""" +pdfsquaredL2norm + +function pdfsquaredL2norm(d::Beta) + α, β = params(d) + z = beta(2 * α - 1, 2 * β - 1) / beta(α, β) ^ 2 + # L2 norm of the pdf converges only for α > 0.5 and β > 0.5 + return α > 0.5 && β > 0.5 ? z : oftype(z, Inf) +end + +pdfsquaredL2norm(d::Cauchy) = inv2π / d.σ + +function pdfsquaredL2norm(d::Chi) + ν = d.ν + z = (2 ^ (1 - ν) * gamma((2 * ν - 1) / 2)) / gamma(ν / 2) ^ 2 + # L2 norm of the pdf converges only for ν > 0.5 + return d.ν > 0.5 ? z : oftype(z, Inf) +end + +function pdfsquaredL2norm(d::Chisq) + ν = d.ν + z = gamma(d.ν - 1) / (gamma(d.ν / 2) ^ 2 * 2 ^ d.ν) + # L2 norm of the pdf converges only for ν > 1 + return ν > 1 ? z : oftype(z, Inf) +end + +pdfsquaredL2norm(d::Exponential) = 1 / (2 * d.θ) + +function pdfsquaredL2norm(d::Gamma) + α, θ = params(d) + z = (2^(1 - 2 * α) * gamma(2 * α - 1)) / (gamma(α) ^ 2 * θ) + # L2 norm of the pdf converges only for α > 0.5 + return α > 0.5 ? z : oftype(z, Inf) +end + +pdfsquaredL2norm(d::Logistic) = 1 / (6 * d.θ) + +pdfsquaredL2norm(d::Normal) = inv(sqrt4π * d.σ) + +pdfsquaredL2norm(d::Uniform) = 1 / (d.b - d.a) diff --git a/test/pdfnorm.jl b/test/pdfnorm.jl new file mode 100644 index 000000000..bf9334910 --- /dev/null +++ b/test/pdfnorm.jl @@ -0,0 +1,63 @@ +using Test, Distributions, SpecialFunctions + +@testset "pdf L2 norm" begin + # Test error on a non implemented norm. + @test_throws MethodError pdfsquaredL2norm(Gumbel()) + + @testset "Beta" begin + @test pdfsquaredL2norm(Beta(1, 1)) ≈ 1 + @test pdfsquaredL2norm(Beta(2, 2)) ≈ 6 / 5 + @test pdfsquaredL2norm(Beta(0.25, 1)) ≈ Inf + @test pdfsquaredL2norm(Beta(1, 0.25)) ≈ Inf + end + + @testset "Cauchy" begin + @test pdfsquaredL2norm(Cauchy(0, 1)) ≈ 1 / (2 * π) + @test pdfsquaredL2norm(Cauchy(0, 2)) ≈ 1 / (4 * π) + # The norm doesn't depend on the mean + @test pdfsquaredL2norm(Cauchy(100, 1)) == pdfsquaredL2norm(Cauchy(-100, 1)) == pdfsquaredL2norm(Cauchy(0, 1)) + end + + @testset "Chi" begin + @test pdfsquaredL2norm(Chi(2)) ≈ gamma(3 / 2) / 2 + @test pdfsquaredL2norm(Chi(0.25)) ≈ Inf + end + + @testset "Chisq" begin + @test pdfsquaredL2norm(Chisq(2)) ≈ 1 / 4 + @test pdfsquaredL2norm(Chisq(1)) ≈ Inf + end + + @testset "Exponential" begin + @test pdfsquaredL2norm(Exponential(1)) ≈ 1 / 2 + @test pdfsquaredL2norm(Exponential(2)) ≈ 1 / 4 + end + + @testset "Gamma" begin + @test pdfsquaredL2norm(Gamma(1, 1)) ≈ 1 / 2 + @test pdfsquaredL2norm(Gamma(1, 2)) ≈ 1 / 4 + @test pdfsquaredL2norm(Gamma(2, 2)) ≈ 1 / 8 + @test pdfsquaredL2norm(Gamma(1, 0.25)) ≈ 2 + @test pdfsquaredL2norm(Gamma(0.5, 1)) ≈ Inf + end + + @testset "Logistic" begin + @test pdfsquaredL2norm(Logistic(0, 1)) ≈ 1 / 6 + @test pdfsquaredL2norm(Logistic(0, 2)) ≈ 1 / 12 + # The norm doesn't depend on the mean + @test pdfsquaredL2norm(Logistic(100, 1)) == pdfsquaredL2norm(Logistic(-100, 1)) == pdfsquaredL2norm(Logistic(0, 1)) + end + + @testset "Normal" begin + @test pdfsquaredL2norm(Normal(0, 1)) ≈ 1 / (2 * sqrt(π)) + @test pdfsquaredL2norm(Normal(0, 2)) ≈ 1 / (4 * sqrt(π)) + @test pdfsquaredL2norm(Normal(1, 0)) ≈ Inf + # The norm doesn't depend on the mean + @test pdfsquaredL2norm(Normal(100, 1)) == pdfsquaredL2norm(Normal(-100, 1)) == pdfsquaredL2norm(Normal(0, 1)) + end + + @testset "Uniform" begin + @test pdfsquaredL2norm(Uniform(-1, 1)) ≈ 1 / 2 + @test pdfsquaredL2norm(Uniform(1, 2)) ≈ 1 + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 2f3cc3f68..e2a83166a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -61,6 +61,7 @@ const tests = [ "skewnormal", "chi", "gumbel", + "pdfnorm", ] printstyled("Running tests:\n", color=:blue)