Add function to compute L2 norm of distributions (#1321)

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

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"

docs/src/univariate.md

@@ -63,6 +63,7 @@ entropy(::UnivariateDistribution, ::Bool)
 entropy(::UnivariateDistribution, ::Real)
 mgf(::UnivariateDistribution, ::Any)
 cf(::UnivariateDistribution, ::Any)
+pdfsquaredL2norm
 ```
 
 ### Probability Evaluation

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")

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)

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

test/runtests.jl

@@ -61,6 +61,7 @@ const tests = [
     "skewnormal",
     "chi",
     "gumbel",
+    "pdfnorm",
 ]
 
 printstyled("Running tests:\n", color=:blue)