Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add squared L2 norms of some discrete distributions #1340

Merged
merged 3 commits into from
Jun 9, 2021
Merged

Add squared L2 norms of some discrete distributions #1340

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d181533
Add squared L2 norms of some discrete distributions
giordano Jun 7, 2021
7641746
Compute L2 norms of PDFs numerically in tests
giordano Jun 8, 2021
91a4ffa
Update src/pdfnorm.jl
giordano Jun 8, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
15 changes: 15 additions & 0 deletions src/pdfnorm.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
giordano marked this conversation as resolved.
Show resolved Hide resolved

pdfsquaredL2norm(d::Cauchy) = inv2π / d.σ

function pdfsquaredL2norm(d::Chi)
Expand All @@ -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)
Expand All @@ -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>).
Comment on lines +58 to +62
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Computing this took me more time than I'd like to admit. If you have the definition of the modified Bessel function of the first kind, then the identity is easy to see, but I had to dig it up quite a bit to get there. I don't know whether the reference I dropped here is more useful or confusing, I can remove it if you consider it unnecessary

pdfsquaredL2norm(d::Poisson) = besseli(0, 2 * d.λ) * exp(-2 * d.λ)

pdfsquaredL2norm(d::Uniform) = 1 / (d.b - d.a)
92 changes: 74 additions & 18 deletions test/pdfnorm.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,63 +1,119 @@
using Test, Distributions, SpecialFunctions
using QuadGK

# `numeric_norm` is a helper function to compute numerically the squared L2
# norms of the distributions. These methods aren't very robust because can't
# deal with divergent norms, or discrete distributions with infinite support.
numeric_norm(d::ContinuousUnivariateDistribution) =
quadgk(x -> pdf(d, x) ^ 2, support(d).lb, support(d).ub)[1]

function numeric_norm(d::DiscreteUnivariateDistribution)
# When the distribution has infinite support, sum up to an arbitrary large
# value.
upper = isfinite(maximum(d)) ? round(Int, maximum(d)) : 100
return sum(pdf(d, k) ^ 2 for k in round(Int, minimum(d)):upper)
end

@testset "pdf L2 norm" begin
# Test error on a non implemented norm.
@test_throws MethodError pdfsquaredL2norm(Gumbel())

@testset "Bernoulli" begin
for d in (Bernoulli(0.5), Bernoulli(0), Bernoulli(0.25))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
# The norm is the same for complementary probabilities
@test pdfsquaredL2norm(Bernoulli(0)) == pdfsquaredL2norm(Bernoulli(1))
@test pdfsquaredL2norm(Bernoulli(0.25)) == pdfsquaredL2norm(Bernoulli(0.75))
end

@testset "Beta" begin
@test pdfsquaredL2norm(Beta(1, 1)) ≈ 1
@test pdfsquaredL2norm(Beta(2, 2)) ≈ 6 / 5
for d in (Beta(1, 1), Beta(2, 2))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
@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)
d = Categorical(collect(1 / n for _ in 1:n))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
for d in (Categorical([0.25, 0.75]), Categorical([1 / 6, 1 / 3, 1 / 2]))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
end

@testset "Cauchy" begin
@test pdfsquaredL2norm(Cauchy(0, 1)) ≈ 1 / (2 * π)
@test pdfsquaredL2norm(Cauchy(0, 2)) ≈ 1 / (4 * π)
for d in (Cauchy(0, 1), Cauchy(0, 2))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
# 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(2)) ≈ numeric_norm(Chi(2))
@test pdfsquaredL2norm(Chi(0.25)) ≈ Inf
end

@testset "Chisq" begin
@test pdfsquaredL2norm(Chisq(2)) ≈ 1 / 4
@test pdfsquaredL2norm(Chisq(2)) ≈ numeric_norm(Chisq(2))
@test pdfsquaredL2norm(Chisq(1)) ≈ Inf
end

@testset "DiscreteUniform" begin
for d in (DiscreteUniform(-1, 1), DiscreteUniform(1, 2))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
end

@testset "Exponential" begin
@test pdfsquaredL2norm(Exponential(1)) ≈ 1 / 2
@test pdfsquaredL2norm(Exponential(2)) ≈ 1 / 4
for d in (Exponential(1), Exponential(2))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
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
for d in (Gamma(1, 1), Gamma(1, 2), Gamma(2, 2), Gamma(1, 0.25))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
@test pdfsquaredL2norm(Gamma(0.5, 1)) ≈ Inf
end

@testset "Geometric" begin
for d in (Geometric(0.20), Geometric(0.25), Geometric(0.50), Geometric(0.75), Geometric(0.80))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
end

@testset "Logistic" begin
@test pdfsquaredL2norm(Logistic(0, 1)) ≈ 1 / 6
@test pdfsquaredL2norm(Logistic(0, 2)) ≈ 1 / 12
for d in (Logistic(0, 1), Logistic(0, 2))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
# 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(π))
for d in (Normal(0, 1), Normal(0, 2))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
@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 "Poisson" begin
for d in (Poisson(0), Poisson(1), Poisson(pi))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
end

@testset "Uniform" begin
@test pdfsquaredL2norm(Uniform(-1, 1)) ≈ 1 / 2
@test pdfsquaredL2norm(Uniform(1, 2)) ≈ 1
for d in (Uniform(-1, 1), Uniform(1, 2))
@test pdfsquaredL2norm(d) ≈ numeric_norm(d)
end
end
end