Add squared L2 norms of some discrete distributions

Project.toml

@@ -1,7 +1,7 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.2"
+version = "0.25.3"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

src/pdfnorm.jl

@@ -11,13 +11,17 @@ where ``S`` is the support of ``f(x)``.
 """
 pdfsquaredL2norm
 
+pdfsquaredL2norm(d::Bernoulli) = @evalpoly d.p 1 -2 2
+
 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::Categorical) = dot(probs(d), probs(d))
+
 pdfsquaredL2norm(d::Cauchy) = inv2π / d.σ
 
 function pdfsquaredL2norm(d::Chi)
@@ -34,6 +38,8 @@ function pdfsquaredL2norm(d::Chisq)
     return ν > 1 ? z : oftype(z, Inf)
 end
 
+pdfsquaredL2norm(d::DiscreteUniform) = 1 / (d.b - d.a + 1)
+
 pdfsquaredL2norm(d::Exponential) = 1 / (2 * d.θ)
 
 function pdfsquaredL2norm(d::Gamma)
@@ -43,8 +49,17 @@ function pdfsquaredL2norm(d::Gamma)
     return α > 0.5 ? z : oftype(z, Inf)
 end
 
+pdfsquaredL2norm(d::Geometric) = d.p ^ 2 / (2 * d.p - d.p ^ 2)
+
 pdfsquaredL2norm(d::Logistic) = 1 / (6 * d.θ)
 
 pdfsquaredL2norm(d::Normal) = inv(sqrt4π * d.σ)
 
+# The identity is obvious if you look at the definition of the modified Bessel function of
+# first kind I_0.  Starting from the L2-norm of the Poisson distribution this can be proven
+# by observing this is the Laguerre exponential l−e of x², which is related to said Bessel
+# function, see <https://doi.org/10.1140/epjst/e2018-00073-1> (preprint:
+# <https://arxiv.org/abs/1707.01135>).
+pdfsquaredL2norm(d::Poisson) = besseli(0, 2 * d.λ) * exp(-2 * d.λ)
+
 pdfsquaredL2norm(d::Uniform) = 1 / (d.b - d.a)

test/pdfnorm.jl

@@ -4,13 +4,28 @@ using Test, Distributions, SpecialFunctions
     # Test error on a non implemented norm.
     @test_throws MethodError pdfsquaredL2norm(Gumbel())
 
+    @testset "Bernoulli" begin
+        @test pdfsquaredL2norm(Bernoulli(0.5)) ≈ 0.5
+        # The norm is the same for complementary probabilities
+        @test pdfsquaredL2norm(Bernoulli(0)) == pdfsquaredL2norm(Bernoulli(1)) ≈ 1
+        @test pdfsquaredL2norm(Bernoulli(0.25)) == pdfsquaredL2norm(Bernoulli(0.75)) ≈ 0.625
+    end
+
     @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 "Categorical" begin
+        for n in (1, 2, 5, 10)
+            @test pdfsquaredL2norm(Categorical(collect(1 / n for _ in 1:n))) ≈ 1 / n
+        end
+        @test pdfsquaredL2norm(Categorical([0.25, 0.75])) ≈ 0.625
+        @test pdfsquaredL2norm(Categorical([1 / 6, 1 / 3, 1 / 2])) ≈ 7 / 18
+    end
+
     @testset "Cauchy" begin
         @test pdfsquaredL2norm(Cauchy(0, 1)) ≈ 1 / (2 * π)
         @test pdfsquaredL2norm(Cauchy(0, 2)) ≈ 1 / (4 * π)
@@ -28,6 +43,11 @@ using Test, Distributions, SpecialFunctions
         @test pdfsquaredL2norm(Chisq(1)) ≈ Inf
     end
 
+    @testset "DiscreteUniform" begin
+        @test pdfsquaredL2norm(DiscreteUniform(-1, 1)) ≈ 1 / 3
+        @test pdfsquaredL2norm(DiscreteUniform(1, 2)) ≈ 1 / 2
+    end
+
     @testset "Exponential" begin
         @test pdfsquaredL2norm(Exponential(1)) ≈ 1 / 2
         @test pdfsquaredL2norm(Exponential(2)) ≈ 1 / 4
@@ -41,6 +61,14 @@ using Test, Distributions, SpecialFunctions
         @test pdfsquaredL2norm(Gamma(0.5, 1)) ≈ Inf
     end
 
+    @testset "Geometric" begin
+        @test pdfsquaredL2norm(Geometric(0.20)) ≈ 1 / 9
+        @test pdfsquaredL2norm(Geometric(0.25)) ≈ 1 / 7
+        @test pdfsquaredL2norm(Geometric(0.50)) ≈ 1 / 3
+        @test pdfsquaredL2norm(Geometric(0.75)) ≈ 3 / 5
+        @test pdfsquaredL2norm(Geometric(0.80)) ≈ 2 / 3
+    end
+
     @testset "Logistic" begin
         @test pdfsquaredL2norm(Logistic(0, 1)) ≈ 1 / 6
         @test pdfsquaredL2norm(Logistic(0, 2)) ≈ 1 / 12
@@ -56,6 +84,12 @@ using Test, Distributions, SpecialFunctions
         @test pdfsquaredL2norm(Normal(100, 1)) == pdfsquaredL2norm(Normal(-100, 1)) == pdfsquaredL2norm(Normal(0, 1))
     end
 
+    @testset "Poisson" begin
+        @test pdfsquaredL2norm(Poisson(0)) ≈ 1
+        @test pdfsquaredL2norm(Poisson(1)) ≈ besseli(0, 2) * exp(-2)
+        @test pdfsquaredL2norm(Poisson(pi)) ≈ besseli(0, 2 * pi) * exp(-2)
+    end
+
     @testset "Uniform" begin
         @test pdfsquaredL2norm(Uniform(-1, 1)) ≈ 1 / 2
         @test pdfsquaredL2norm(Uniform(1, 2)) ≈ 1