From 932f5165bfc6997118b797336f6e6941c7ef86df Mon Sep 17 00:00:00 2001 From: st-- Date: Mon, 28 Mar 2022 14:49:55 +0200 Subject: [PATCH 01/44] less code duplication (#72) Co-authored-by: David Widmann --- README.md | 8 ++++++-- src/likelihoods/bernoulli.jl | 2 +- src/likelihoods/categorical.jl | 4 +--- src/likelihoods/exponential.jl | 4 +--- src/likelihoods/gamma.jl | 2 +- src/likelihoods/negativebinomial.jl | 4 +--- src/likelihoods/poisson.jl | 2 +- 7 files changed, 12 insertions(+), 14 deletions(-) diff --git a/README.md b/README.md index ae9027c..d0e86a3 100644 --- a/README.md +++ b/README.md @@ -7,5 +7,9 @@ [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) [![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac) -GPLikelihoods.jl provides a collection of likelihoods to be used as building blocks for defining non-Gaussian problems. -It is intended to be mainly consumed by other packages within the JuliaGaussianProcesses ecosystem (such as [ApproximateGPs.jl](https://github.com/JuliaGaussianProcesses/ApproximateGPs.jl)) for approximate Gaussian process regression, classification, event counting, etc. +GPLikelihoods.jl provides a collection of likelihoods to be used as building +blocks for defining non-Gaussian problems. It is intended to be mainly +consumed by other packages within the JuliaGaussianProcesses ecosystem (such as +[ApproximateGPs.jl](https://github.com/JuliaGaussianProcesses/ApproximateGPs.jl)) +for approximate Gaussian process regression, classification, event counting, +etc. diff --git a/src/likelihoods/bernoulli.jl b/src/likelihoods/bernoulli.jl index 3866cb5..d14878e 100644 --- a/src/likelihoods/bernoulli.jl +++ b/src/likelihoods/bernoulli.jl @@ -18,4 +18,4 @@ BernoulliLikelihood(l=logistic) = BernoulliLikelihood(link(l)) (l::BernoulliLikelihood)(f::Real) = Bernoulli(l.invlink(f)) -(l::BernoulliLikelihood)(fs::AbstractVector{<:Real}) = Product(Bernoulli.(l.invlink.(fs))) +(l::BernoulliLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) diff --git a/src/likelihoods/categorical.jl b/src/likelihoods/categorical.jl index 10c5c36..6caedbf 100644 --- a/src/likelihoods/categorical.jl +++ b/src/likelihoods/categorical.jl @@ -36,6 +36,4 @@ function (l::CategoricalLikelihood)(f::AbstractVector{<:Real}) return Categorical(l.invlink(f)) end -function (l::CategoricalLikelihood)(fs::AbstractVector) - return Product(Categorical.(l.invlink.(fs))) -end +(l::CategoricalLikelihood)(fs::AbstractVector) = Product(map(l, fs)) diff --git a/src/likelihoods/exponential.jl b/src/likelihoods/exponential.jl index fdfacb8..ee1483e 100644 --- a/src/likelihoods/exponential.jl +++ b/src/likelihoods/exponential.jl @@ -15,6 +15,4 @@ ExponentialLikelihood(l=exp) = ExponentialLikelihood(link(l)) (l::ExponentialLikelihood)(f::Real) = Exponential(l.invlink(f)) -function (l::ExponentialLikelihood)(fs::AbstractVector{<:Real}) - return Product(Exponential.(l.invlink.(fs))) -end +(l::ExponentialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) diff --git a/src/likelihoods/gamma.jl b/src/likelihoods/gamma.jl index bd13dea..512d56a 100644 --- a/src/likelihoods/gamma.jl +++ b/src/likelihoods/gamma.jl @@ -20,4 +20,4 @@ GammaLikelihood(α::Real=1.0, l=exp) = GammaLikelihood(α, link(l)) (l::GammaLikelihood)(f::Real) = Gamma(l.α, l.invlink(f)) -(l::GammaLikelihood)(fs::AbstractVector{<:Real}) = Product(Gamma.(l.α, l.invlink.(fs))) +(l::GammaLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index a48b7f2..0b6477d 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -33,6 +33,4 @@ end (l::NegativeBinomialLikelihood)(f::Real) = NegativeBinomial(l.successes, l.invlink(f)) -function (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) - return Product(NegativeBinomial.(l.successes, l.invlink.(fs))) -end +(l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) diff --git a/src/likelihoods/poisson.jl b/src/likelihoods/poisson.jl index c8f7e28..f8b7b50 100644 --- a/src/likelihoods/poisson.jl +++ b/src/likelihoods/poisson.jl @@ -18,4 +18,4 @@ PoissonLikelihood(l=exp) = PoissonLikelihood(link(l)) (l::PoissonLikelihood)(f::Real) = Poisson(l.invlink(f)) -(l::PoissonLikelihood)(fs::AbstractVector{<:Real}) = Product(Poisson.(l.invlink.(fs))) +(l::PoissonLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) From 18ae4362eda2f2218c6fcdc84051dea306505247 Mon Sep 17 00:00:00 2001 From: st-- Date: Mon, 28 Mar 2022 16:42:45 +0200 Subject: [PATCH 02/44] move expected_loglik here from ApproximateGPs.jl (#70) Moves `expected_loglik` from ApproximateGPs.jl to GPLikelihoods.jlApproximateGPs. Changes from the original code in ApproximateGPs.jl: - rename to `expected_loglikelihood` to be more explicit - changed signature: argument order is now `quadrature`, `lik`, `q_f`, `y` - rename `_default_quadrature` to `default_expectation_method` - rename expectation structs to `DefaultExpectationMethod`, `AnalyticExpectation`, `GaussHermiteExpectation`, `MonteCarloExpectation`; they are not exported from here - drop abstract base type (was not used for anything) - drop default argument for GaussHermite / MonteCarlo quadrature points (moved into `DefaultExpectationMethod`) - GaussHermiteExpectation now stores nodes&weights within struct, resolves #71 - move analytic definitions into files for each likelihood implementation Note: bumps julia compat to 1.6. Co-authored-by: willtebbutt --- .github/workflows/CI.yml | 2 +- Project.toml | 12 +++- src/GPLikelihoods.jl | 2 + src/expectations.jl | 117 +++++++++++++++++++++++++++++++++ src/likelihoods/exponential.jl | 12 ++++ src/likelihoods/gamma.jl | 15 +++++ src/likelihoods/gaussian.jl | 13 ++++ src/likelihoods/poisson.jl | 12 ++++ test/Project.toml | 1 + test/expectations.jl | 100 ++++++++++++++++++++++++++++ test/runtests.jl | 3 + 11 files changed, 286 insertions(+), 3 deletions(-) create mode 100644 src/expectations.jl create mode 100644 test/expectations.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 2154e59..99717df 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -22,7 +22,7 @@ jobs: matrix: version: - '1' - - '1.3' + - '1.6' - 'nightly' os: - ubuntu-latest diff --git a/Project.toml b/Project.toml index 7e8c25e..8ee087a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,20 +1,28 @@ name = "GPLikelihoods" uuid = "6031954c-0455-49d7-b3b9-3e1c99afaf40" authors = ["JuliaGaussianProcesses Team"] -version = "0.3.1" +version = "0.4.0" [deps] +ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196" InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" +IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] +ChainRulesCore = "1.7" Distributions = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25" +FastGaussQuadrature = "0.4" Functors = "0.1, 0.2" InverseFunctions = "0.1.2" +IrrationalConstants = "0.1" +SpecialFunctions = "1, 2" StatsFuns = "0.9.13" -julia = "1.3" +julia = "1.6" diff --git a/src/GPLikelihoods.jl b/src/GPLikelihoods.jl index 3eb40a4..33973cb 100644 --- a/src/GPLikelihoods.jl +++ b/src/GPLikelihoods.jl @@ -34,6 +34,8 @@ include("links.jl") # Likelihoods abstract type AbstractLikelihood end + +include("expectations.jl") include("likelihoods/bernoulli.jl") include("likelihoods/categorical.jl") include("likelihoods/gaussian.jl") diff --git a/src/expectations.jl b/src/expectations.jl new file mode 100644 index 0000000..8855cfb --- /dev/null +++ b/src/expectations.jl @@ -0,0 +1,117 @@ +using FastGaussQuadrature: gausshermite +using SpecialFunctions: loggamma +using ChainRulesCore: ChainRulesCore +using IrrationalConstants: sqrt2, invsqrtπ + +struct DefaultExpectationMethod end + +struct AnalyticExpectation end + +struct GaussHermiteExpectation + xs::Vector{Float64} + ws::Vector{Float64} +end +GaussHermiteExpectation(n::Integer) = GaussHermiteExpectation(gausshermite(n)...) + +ChainRulesCore.@non_differentiable gausshermite(n) + +struct MonteCarloExpectation + n_samples::Int +end + +default_expectation_method(_) = GaussHermiteExpectation(20) + +""" + expected_loglikelihood( + quadrature, + lik, + q_f::AbstractVector{<:Normal}, + y::AbstractVector, + ) + +This function computes the expected log likelihood: + +```math + ∫ q(f) log p(y | f) df +``` +where `p(y | f)` is the process likelihood. This is described by `lik`, which should be a +callable that takes `f` as input and returns a Distribution over `y` that supports +`loglikelihood(lik(f), y)`. + +`q(f)` is an approximation to the latent function values `f` given by: +```math + q(f) = ∫ p(f | u) q(u) du +``` +where `q(u)` is the variational distribution over inducing points (see +[`elbo`](@ref)). The marginal distributions of `q(f)` are given by `q_f`. + +`quadrature` determines which method is used to calculate the expected log +likelihood - see [`elbo`](@ref) for more details. + +# Extended help + +`q(f)` is assumed to be an `MvNormal` distribution and `p(y | f)` is assumed to +have independent marginals such that only the marginals of `q(f)` are required. +""" +expected_loglikelihood(quadrature, lik, q_f, y) + +""" + expected_loglikelihood(::DefaultExpectationMethod, lik, q_f::AbstractVector{<:Normal}, y::AbstractVector) + +The expected log likelihood, using the default quadrature method for the given likelihood. +(The default quadrature method is defined by `default_expectation_method(lik)`, and should +be the closed form solution if it exists, but otherwise defaults to Gauss-Hermite +quadrature.) +""" +function expected_loglikelihood( + ::DefaultExpectationMethod, lik, q_f::AbstractVector{<:Normal}, y::AbstractVector +) + quadrature = default_expectation_method(lik) + return expected_loglikelihood(quadrature, lik, q_f, y) +end + +function expected_loglikelihood( + mc::MonteCarloExpectation, lik, q_f::AbstractVector{<:Normal}, y::AbstractVector +) + # take `n_samples` reparameterised samples + f_μ = mean.(q_f) + fs = f_μ .+ std.(q_f) .* randn(eltype(f_μ), length(q_f), mc.n_samples) + lls = loglikelihood.(lik.(fs), y) + return sum(lls) / mc.n_samples +end + +# Compute the expected_loglikelihood over a collection of observations and marginal distributions +function expected_loglikelihood( + gh::GaussHermiteExpectation, lik, q_f::AbstractVector{<:Normal}, y::AbstractVector +) + # Compute the expectation via Gauss-Hermite quadrature + # using a reparameterisation by change of variable + # (see e.g. en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature) + return sum(Broadcast.instantiate( + Broadcast.broadcasted(y, q_f) do yᵢ, q_fᵢ # Loop over every pair + # of marginal distribution q(fᵢ) and observation yᵢ + expected_loglikelihood(gh, lik, q_fᵢ, yᵢ) + end, + )) +end + +# Compute the expected_loglikelihood for one observation and a marginal distributions +function expected_loglikelihood(gh::GaussHermiteExpectation, lik, q_f::Normal, y) + μ = mean(q_f) + σ̃ = sqrt2 * std(q_f) + return invsqrtπ * sum(Broadcast.instantiate( + Broadcast.broadcasted(gh.xs, gh.ws) do x, w # Loop over every + # pair of Gauss-Hermite point x with weight w + f = σ̃ * x + μ + loglikelihood(lik(f), y) * w + end, + )) +end + +function expected_loglikelihood( + ::AnalyticExpectation, lik, q_f::AbstractVector{<:Normal}, y::AbstractVector +) + return error( + "No analytic solution exists for $(typeof(lik)). Use `DefaultExpectationMethod`, `GaussHermiteExpectation` or `MonteCarloExpectation` instead.", + ) +end diff --git a/src/likelihoods/exponential.jl b/src/likelihoods/exponential.jl index ee1483e..1b7dc60 100644 --- a/src/likelihoods/exponential.jl +++ b/src/likelihoods/exponential.jl @@ -16,3 +16,15 @@ ExponentialLikelihood(l=exp) = ExponentialLikelihood(link(l)) (l::ExponentialLikelihood)(f::Real) = Exponential(l.invlink(f)) (l::ExponentialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) + +function expected_loglikelihood( + ::AnalyticExpectation, + ::ExponentialLikelihood{ExpLink}, + q_f::AbstractVector{<:Normal}, + y::AbstractVector{<:Real}, +) + f_μ = mean.(q_f) + return sum(-f_μ - y .* exp.((var.(q_f) / 2) .- f_μ)) +end + +default_expectation_method(::ExponentialLikelihood{ExpLink}) = AnalyticExpectation() diff --git a/src/likelihoods/gamma.jl b/src/likelihoods/gamma.jl index 512d56a..dcd0202 100644 --- a/src/likelihoods/gamma.jl +++ b/src/likelihoods/gamma.jl @@ -21,3 +21,18 @@ GammaLikelihood(α::Real=1.0, l=exp) = GammaLikelihood(α, link(l)) (l::GammaLikelihood)(f::Real) = Gamma(l.α, l.invlink(f)) (l::GammaLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) + +function expected_loglikelihood( + ::AnalyticExpectation, + lik::GammaLikelihood{ExpLink}, + q_f::AbstractVector{<:Normal}, + y::AbstractVector{<:Real}, +) + f_μ = mean.(q_f) + return sum( + (lik.α - 1) * log.(y) .- y .* exp.((var.(q_f) / 2) .- f_μ) .- lik.α * f_μ .- + loggamma(lik.α), + ) +end + +default_expectation_method(::GammaLikelihood{ExpLink}) = AnalyticExpectation() diff --git a/src/likelihoods/gaussian.jl b/src/likelihoods/gaussian.jl index a20fa26..ca5101b 100644 --- a/src/likelihoods/gaussian.jl +++ b/src/likelihoods/gaussian.jl @@ -23,6 +23,19 @@ GaussianLikelihood(σ²::Real) = GaussianLikelihood([σ²]) (l::GaussianLikelihood)(fs::AbstractVector{<:Real}) = MvNormal(fs, first(l.σ²) * I) +function expected_loglikelihood( + ::AnalyticExpectation, + lik::GaussianLikelihood, + q_f::AbstractVector{<:Normal}, + y::AbstractVector{<:Real}, +) + return sum( + -0.5 * (log(2π) .+ log.(lik.σ²) .+ ((y .- mean.(q_f)) .^ 2 .+ var.(q_f)) / lik.σ²) + ) +end + +default_expectation_method(::GaussianLikelihood) = AnalyticExpectation() + """ HeteroscedasticGaussianLikelihood(l=exp) diff --git a/src/likelihoods/poisson.jl b/src/likelihoods/poisson.jl index f8b7b50..9e31d6e 100644 --- a/src/likelihoods/poisson.jl +++ b/src/likelihoods/poisson.jl @@ -19,3 +19,15 @@ PoissonLikelihood(l=exp) = PoissonLikelihood(link(l)) (l::PoissonLikelihood)(f::Real) = Poisson(l.invlink(f)) (l::PoissonLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) + +function expected_loglikelihood( + ::AnalyticExpectation, + ::PoissonLikelihood{ExpLink}, + q_f::AbstractVector{<:Normal}, + y::AbstractVector{<:Real}, +) + f_μ = mean.(q_f) + return sum((y .* f_μ) - exp.(f_μ .+ (var.(q_f) / 2)) - loggamma.(y .+ 1)) +end + +default_expectation_method(::PoissonLikelihood{ExpLink}) = AnalyticExpectation() diff --git a/test/Project.toml b/test/Project.toml index 259ad91..b88386d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,7 @@ Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] Distributions = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25" diff --git a/test/expectations.jl b/test/expectations.jl new file mode 100644 index 0000000..65fe96e --- /dev/null +++ b/test/expectations.jl @@ -0,0 +1,100 @@ +@testset "expectations" begin + # Test that the various methods of computing expectations return the same + # result. + rng = MersenneTwister(123456) + q_f = Normal.(zeros(10), ones(10)) + + likelihoods_to_test = [ + ExponentialLikelihood(), + GammaLikelihood(), + PoissonLikelihood(), + GaussianLikelihood(), + ] + + @testset "testing all analytic implementations" begin + # Test that we're not missing any analytic implementation in `likelihoods_to_test`! + implementation_types = [ + (; quadrature=m.sig.types[2], lik=m.sig.types[3]) for + m in methods(GPLikelihoods.expected_loglikelihood) + ] + analytic_likelihoods = [ + m.lik for m in implementation_types if + m.quadrature == GPLikelihoods.AnalyticExpectation && m.lik != Any + ] + for lik_type in analytic_likelihoods + lik_type_instances = filter(lik -> isa(lik, lik_type), likelihoods_to_test) + @test !isempty(lik_type_instances) + lik = first(lik_type_instances) + @test GPLikelihoods.default_expectation_method(lik) isa + GPLikelihoods.AnalyticExpectation + end + end + + @testset "$(nameof(typeof(lik)))" for lik in likelihoods_to_test + methods = [ + GaussHermiteExpectation(100), + MonteCarloExpectation(1e7), + GPLikelihoods.DefaultExpectationMethod(), + ] + def = GPLikelihoods.default_expectation_method(lik) + if def isa GPLikelihoods.AnalyticExpectation + push!(methods, def) + end + y = rand.(rng, lik.(zeros(10))) + + results = map(m -> GPLikelihoods.expected_loglikelihood(m, lik, q_f, y), methods) + @test all(x -> isapprox(x, results[end]; atol=1e-6, rtol=1e-3), results) + end + + @test GPLikelihoods.expected_loglikelihood( + MonteCarloExpectation(1), GaussianLikelihood(), q_f, zeros(10) + ) isa Real + @test GPLikelihoods.expected_loglikelihood( + GaussHermiteExpectation(1), GaussianLikelihood(), q_f, zeros(10) + ) isa Real + @test GPLikelihoods.default_expectation_method(θ -> Normal(0, θ)) isa + GaussHermiteExpectation + + # see https://github.com/JuliaGaussianProcesses/ApproximateGPs.jl/issues/82 + @testset "testing Zygote compatibility with GaussHermiteExpectation" begin + N = 10 + gh = GaussHermiteExpectation(12) + μs = randn(rng, N) + σs = rand(rng, N) + + # Test differentiation with variational parameters + for lik in likelihoods_to_test + y = rand.(rng, lik.(rand.(Normal.(μs, σs)))) + gμ, glogσ = Zygote.gradient(μs, log.(σs)) do μ, logσ + GPLikelihoods.expected_loglikelihood(gh, lik, Normal.(μ, exp.(logσ)), y) + end + @test all(isfinite, gμ) + @test all(isfinite, glogσ) + end + + # Test differentiation with likelihood parameters + # Test GaussianLikelihood parameter + σ = 1.0 + y = randn(rng, N) + glogσ = only( + Zygote.gradient(log(σ)) do x + GPLikelihoods.expected_loglikelihood( + gh, GaussianLikelihood(exp(x)), Normal.(μs, σs), y + ) + end, + ) + @test isfinite(glogσ) + + # Test GammaLikelihood parameter + α = 2.0 + y = rand.(rng, Gamma.(α, rand(N))) + glogα = only( + Zygote.gradient(log(α)) do x + GPLikelihoods.expected_loglikelihood( + gh, GammaLikelihood(exp(x)), Normal.(μs, σs), y + ) + end, + ) + @test isfinite(glogα) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index d8a582e..1552a99 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,10 +1,12 @@ using GPLikelihoods +using GPLikelihoods: GaussHermiteExpectation, MonteCarloExpectation using GPLikelihoods.TestInterface: test_interface using Test using Random using Functors using Distributions using StatsFuns +using Zygote @testset "GPLikelihoods.jl" begin include("links.jl") @@ -17,4 +19,5 @@ using StatsFuns include("likelihoods/exponential.jl") include("likelihoods/negativebinomial.jl") end + include("expectations.jl") end From 04155ec717160b876f8cba67f0d8e89e7c7f7f43 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Tue, 29 Mar 2022 09:07:10 +0300 Subject: [PATCH 03/44] CompatHelper: add new compat entry for "Zygote" at version "0.6" for package test (#74) Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index b88386d..9df0ed1 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -10,4 +10,5 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" Distributions = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25" Functors = "0.1, 0.2" StatsFuns = "0.9" +Zygote = "0.6" julia = "1.3" From 498af318aecb891c1497e1a4a19c26dc7ebd0c47 Mon Sep 17 00:00:00 2001 From: st-- Date: Tue, 29 Mar 2022 21:40:42 +0200 Subject: [PATCH 04/44] export expected_loglikelihood (#75) --- Project.toml | 2 +- src/GPLikelihoods.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8ee087a..f5b9c6d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GPLikelihoods" uuid = "6031954c-0455-49d7-b3b9-3e1c99afaf40" authors = ["JuliaGaussianProcesses Team"] -version = "0.4.0" +version = "0.4.1" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" diff --git a/src/GPLikelihoods.jl b/src/GPLikelihoods.jl index 33973cb..67e3583 100644 --- a/src/GPLikelihoods.jl +++ b/src/GPLikelihoods.jl @@ -28,6 +28,7 @@ export Link, ProbitLink, NormalCDFLink, SoftMaxLink +export expected_loglikelihood # Links include("links.jl") From 431918053f82329c16333b189ee451a25ed9d5f9 Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Sun, 17 Apr 2022 21:57:45 +0200 Subject: [PATCH 05/44] Implement multiple parametrization for NegativeBinomial --- docs/src/api.md | 9 ++ src/GPLikelihoods.jl | 1 + src/likelihoods/negativebinomial.jl | 125 +++++++++++++++++++++------ test/likelihoods/negativebinomial.jl | 31 ++++--- 4 files changed, 129 insertions(+), 37 deletions(-) diff --git a/docs/src/api.md b/docs/src/api.md index 48ec7ce..1890a4d 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -15,6 +15,15 @@ HeteroscedasticGaussianLikelihood PoissonLikelihood ``` +### Negative Binomial + +```@docs +NegativeBinomialLikelihood +NBParamI +NBParamII +NBParamIII +``` + ## Links ```@docs diff --git a/src/GPLikelihoods.jl b/src/GPLikelihoods.jl index 67e3583..38fce29 100644 --- a/src/GPLikelihoods.jl +++ b/src/GPLikelihoods.jl @@ -15,6 +15,7 @@ export BernoulliLikelihood, ExponentialLikelihood, GammaLikelihood, NegativeBinomialLikelihood +export NBParamI, NBParamII, NBParamIII export Link, ChainLink, BijectiveSimplexLink, diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 0b6477d..7a97c60 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -1,36 +1,111 @@ +abstract type NBParam end + """ - NegativeBinomialLikelihood(l=logistic; successes::Real=1) + NegativeBinomialLikelihood{NBParam}(l=logistic/exp; kwargs...) -Negative binomial likelihood with number of successes `successes`. +There are many possible parametrizations for the Negative Binomial likelihood. +The `NegativeBinomialLikelihood` has a special structure, the first type `NBParam` +defines what parametrization is used. +Each `NBParam` has its own documentation: +- [`NBParamI`](@ref): This is the definition used by Distributions.jl. +- [`NBParamII`](@ref): This is the definition used by Wikipedia. +- [`NBParamIII`](@ref): This is the definition based on the mean. -```math - p(k|successes, f) = \\frac{\\Gamma(k+successes)}{k! \\Gamma(successes)} l(f)^successes (1 - l(f))^k -``` -On calling, this returns a negative binomial distribution with `successes` successes and -probability of success equal to `l(f)`. - -!!! warning "Parameterization" - The parameter `successes` is different from the parameter `r` in the - [Wikipedia definition](http://en.wikipedia.org/wiki/Negative_binomial_distribution), - which denotes the number of failures. - This parametrization is used in order to stay consistent with the parametrization in - [Distributions.jl](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.NegativeBinomial). - To use the Wikipedia definition, set `successes` as the number of "failures" and - change the probability of success from `l(f)` to `1 - l(f)`. - Note that with symmetric functions like the [`LogisticLink`](@ref), this corresponds to - using `l(-f)`. -""" -struct NegativeBinomialLikelihood{Tl<:AbstractLink,T<:Real} <: AbstractLikelihood - successes::T # number of successes parameter +To create a new parametrization, you need to; +- Create a new typee from `struct MyNBParam <: NBParam end` +- Dispatch [`to_dist_params`](@ref): `to_dist_params(l::NegativeBinomialLikelihood{MyNBParam}, f::Real)` +to return the parameters for the [`NegativeBinomial`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial) constructor +- Write a constructor for `NegativeBinomialLikelihood{MyNBParam}` +""" +struct NegativeBinomialLikelihood{NBParam,Tl<:AbstractLink,T} <: AbstractLikelihood + params::T # Likelihood parameters (depends of NBParam) invlink::Tl + function NegativeBinomialLikelihood{Tparam}(params, invlink) where {Tparam} + invlink = link(invlink) + return new{Tparam,typeof(invlink),typeof(params)}(params, invlink) + end end -function NegativeBinomialLikelihood(l=logistic; successes::Real=1) - return NegativeBinomialLikelihood(successes, link(l)) -end +@deprecate NegativeBinomialLikelihood(l=logistic; successes=1) NegativeBinomialLikelihood{ + NBParamI +}( + l; successes +) @functor NegativeBinomialLikelihood -(l::NegativeBinomialLikelihood)(f::Real) = NegativeBinomial(l.successes, l.invlink(f)) +function (l::NegativeBinomialLikelihood)(f::Real) + return NegativeBinomial(to_dist_params(l, f)...) +end (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) + +""" + to_dist_params(l::NegativeBinomialLikelihood{NBParam}, f::Real)->Tuple{Real,Real} + +Take parameters from a given parametrization `NBParam` to the parametrization +used by `Distributions.jl` +""" +to_dist_params + +""" + NBParamI + +Negative Binomial parametrization with `successes` the number of successes and +`l(f)` the probability of `success`. +This corresponds to the definition used by [Distributions.jl](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial). + +```math + p(k|successes, f) = \\frac{\\Gamma(k+successes)}{k! \\Gamma(successes)} l(f)^successes (1 - l(f))^k +``` +""" +struct NBParamI <: NBParam end + +function to_dist_params(l::NegativeBinomialLikelihood{NBParamI}, f::Real) + return l.params.successes, l.invlink(f) +end + +function NegativeBinomialLikelihood{NBParamI}(l=logistic; successes::Real=1) + return NegativeBinomialLikelihood{NBParamI}((; successes), l) +end + +""" + NBParamII + +Negative Binomial parametrization with `failures` the number of failures and +`l(f)` the probability of `success`. +This corresponds to the definition used by [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution). + +```math + p(k|failures, f) = \\frac{\\Gamma(k+failures)}{k! \\Gamma(failure)} l(f)^k (1 - l(f))^{failures} +``` +""" +struct NBParamII <: NBParam end + +function to_dist_params(l::NegativeBinomialLikelihood{NBParamII}, f::Real) + return l.params.failures, 1 - l.invlink(f) +end + +function NegativeBinomialLikelihood{NBParamII}(l=logistic; failures::Real=1) + return NegativeBinomialLikelihood{NBParamII}((; failures), l) +end + +""" + NBParamIII + +Negative Binomial parametrization with mean `μ=l(f)` and number of successes `successes`. +See the definition given in the [Wikipedia article](https://en.wikipedia.org/wiki/Negative_binomial_distribution#Alternative_formulations) +""" +struct NBParamIII <: NBParam end + +function to_dist_params(l::NegativeBinomialLikelihood{NBParamIII}, f::Real) + μ = l.invlink(f) + r = l.params.successes + return r, μ / (μ + r) +end + +function NegativeBinomialLikelihood{NBParamIII}( + l=logistic; successes::Real=1 +) + return NegativeBinomialLikelihood{NBParamIII}((; successes), l) +end diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index e6e0b37..6f559c5 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -1,16 +1,23 @@ @testset "NegativeBinomialLikelihood" begin - for args in ((), (logistic,), (LogisticLink(),)), kwargs in ((), (; successes=1)) - lik = NegativeBinomialLikelihood(args...; kwargs...) - @test lik isa NegativeBinomialLikelihood{LogisticLink,Int} - @test lik.successes == 1 - end + for (nbparam, kwarg) in ((NBParamI, :successes), (NBParamII, :failures), (NBParamIII, :successes)) + @testset "$(nameof(nbparam))" begin + @eval begin + sym_kwarg = $(Meta.quot(kwarg)) + for args in ((logistic,), (LogisticLink(),)), kwargs in ((), (; $(kwarg)=1)) + lik = NegativeBinomialLikelihood{$(nbparam)}(args...; kwargs...) + @test lik isa NegativeBinomialLikelihood{$(nbparam),LogisticLink,<:NamedTuple{<:Any,<:Tuple{Int}}} + end - for args in ((normcdf,), (NormalCDFLink(),)), kwargs in ((; successes=2.0),) - lik = NegativeBinomialLikelihood(args...; kwargs...) - @test lik isa NegativeBinomialLikelihood{NormalCDFLink,Float64} - @test lik.successes == 2 - end + for args in ((normcdf,), (NormalCDFLink(),)), kwargs in ((; $(kwarg)=2.0),) + lik = NegativeBinomialLikelihood{$(nbparam)}(args...; kwargs...) + @test lik isa NegativeBinomialLikelihood{ + $(nbparam),NormalCDFLink,<:NamedTuple{<:Any,<:Tuple{Float64}} + } + end - lik = NegativeBinomialLikelihood() - test_interface(lik, NegativeBinomial; functor_args=(:successes, :invlink)) + lik = NegativeBinomialLikelihood{$(nbparam)}() + test_interface(lik, NegativeBinomial; functor_args=(:params, :invlink)) + end + end + end end From 913592fce5c9cc59b534a33f727cae2f92a57b36 Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Sun, 17 Apr 2022 22:01:06 +0200 Subject: [PATCH 06/44] Change back default for NBParamIII --- src/likelihoods/negativebinomial.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 7a97c60..e05a9b5 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -105,7 +105,7 @@ function to_dist_params(l::NegativeBinomialLikelihood{NBParamIII}, f::Real) end function NegativeBinomialLikelihood{NBParamIII}( - l=logistic; successes::Real=1 + l=exp; successes::Real=1 ) return NegativeBinomialLikelihood{NBParamIII}((; successes), l) end From 6c36d2fdd16f977a81228471ba148788e5be3474 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Sun, 17 Apr 2022 22:08:38 +0200 Subject: [PATCH 07/44] Update test/likelihoods/negativebinomial.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/likelihoods/negativebinomial.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index 6f559c5..4e4e9a7 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -1,5 +1,6 @@ @testset "NegativeBinomialLikelihood" begin - for (nbparam, kwarg) in ((NBParamI, :successes), (NBParamII, :failures), (NBParamIII, :successes)) + for (nbparam, kwarg) in + ((NBParamI, :successes), (NBParamII, :failures), (NBParamIII, :successes)) @testset "$(nameof(nbparam))" begin @eval begin sym_kwarg = $(Meta.quot(kwarg)) From e54833a92231db8cd3a8a6f64cb79c3df1560934 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Sun, 17 Apr 2022 22:08:48 +0200 Subject: [PATCH 08/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/likelihoods/negativebinomial.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index e05a9b5..9b7111e 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -104,8 +104,6 @@ function to_dist_params(l::NegativeBinomialLikelihood{NBParamIII}, f::Real) return r, μ / (μ + r) end -function NegativeBinomialLikelihood{NBParamIII}( - l=exp; successes::Real=1 -) +function NegativeBinomialLikelihood{NBParamIII}(l=exp; successes::Real=1) return NegativeBinomialLikelihood{NBParamIII}((; successes), l) end From b08aec4ebbfaeb5097448dbc84a107884f69e027 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Sun, 17 Apr 2022 22:08:55 +0200 Subject: [PATCH 09/44] Update test/likelihoods/negativebinomial.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/likelihoods/negativebinomial.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index 4e4e9a7..e3da07d 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -6,7 +6,9 @@ sym_kwarg = $(Meta.quot(kwarg)) for args in ((logistic,), (LogisticLink(),)), kwargs in ((), (; $(kwarg)=1)) lik = NegativeBinomialLikelihood{$(nbparam)}(args...; kwargs...) - @test lik isa NegativeBinomialLikelihood{$(nbparam),LogisticLink,<:NamedTuple{<:Any,<:Tuple{Int}}} + @test lik isa NegativeBinomialLikelihood{ + $(nbparam),LogisticLink,<:NamedTuple{<:Any,<:Tuple{Int}} + } end for args in ((normcdf,), (NormalCDFLink(),)), kwargs in ((; $(kwarg)=2.0),) From 72b7129b3d24da769353330082754e5e73df1d3b Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Tue, 19 Apr 2022 09:53:14 +0200 Subject: [PATCH 10/44] Remove to_dist_params --- src/likelihoods/negativebinomial.jl | 30 ++++++++++++----------------- 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 9b7111e..ea47243 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -13,9 +13,10 @@ Each `NBParam` has its own documentation: To create a new parametrization, you need to; - Create a new typee from `struct MyNBParam <: NBParam end` -- Dispatch [`to_dist_params`](@ref): `to_dist_params(l::NegativeBinomialLikelihood{MyNBParam}, f::Real)` -to return the parameters for the [`NegativeBinomial`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial) constructor -- Write a constructor for `NegativeBinomialLikelihood{MyNBParam}` +- Dispatch `(l::NegativeBinomialLikelihood{MyNBParam})(f::Real)`, which return the [`NegativeBinomial`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial) from `Distributions.jl`. +`NegativeBinomial` follows the parametrization [`NBParamI`](@ref), i.e. the first argument is the number of successes +and the second argument is the probability of success. +- Write a constructor for `NegativeBinomialLikelihood{MyNBParam}(link; kwargs...)` """ struct NegativeBinomialLikelihood{NBParam,Tl<:AbstractLink,T} <: AbstractLikelihood params::T # Likelihood parameters (depends of NBParam) @@ -34,19 +35,12 @@ end @functor NegativeBinomialLikelihood -function (l::NegativeBinomialLikelihood)(f::Real) - return NegativeBinomial(to_dist_params(l, f)...) +function (l::NegativeBinomialLikelihood)(::Real) + error("not implemented for type $(typeof(l)). See `NegativeBinomialLikelihood` docs") end (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) -""" - to_dist_params(l::NegativeBinomialLikelihood{NBParam}, f::Real)->Tuple{Real,Real} - -Take parameters from a given parametrization `NBParam` to the parametrization -used by `Distributions.jl` -""" -to_dist_params """ NBParamI @@ -61,8 +55,8 @@ This corresponds to the definition used by [Distributions.jl](https://juliastats """ struct NBParamI <: NBParam end -function to_dist_params(l::NegativeBinomialLikelihood{NBParamI}, f::Real) - return l.params.successes, l.invlink(f) +function (l::NegativeBinomialLikelihood{NBParamI})(f::Real) + return NegativeBinomial(l.params.successes, l.invlink(f)) end function NegativeBinomialLikelihood{NBParamI}(l=logistic; successes::Real=1) @@ -82,8 +76,8 @@ This corresponds to the definition used by [Wikipedia](https://en.wikipedia.org/ """ struct NBParamII <: NBParam end -function to_dist_params(l::NegativeBinomialLikelihood{NBParamII}, f::Real) - return l.params.failures, 1 - l.invlink(f) +function (l::NegativeBinomialLikelihood{NBParamII})(f::Real) + return NegativeBinomial(l.params.failures, 1 - l.invlink(f)) end function NegativeBinomialLikelihood{NBParamII}(l=logistic; failures::Real=1) @@ -98,10 +92,10 @@ See the definition given in the [Wikipedia article](https://en.wikipedia.org/wik """ struct NBParamIII <: NBParam end -function to_dist_params(l::NegativeBinomialLikelihood{NBParamIII}, f::Real) +function (l::NegativeBinomialLikelihood{NBParamIII})(f::Real) μ = l.invlink(f) r = l.params.successes - return r, μ / (μ + r) + return NegativeBinomial(r, μ / (μ + r)) end function NegativeBinomialLikelihood{NBParamIII}(l=exp; successes::Real=1) From a4333908bbab94bf21885ab2c3434e5e66330009 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 19 Apr 2022 10:20:49 +0200 Subject: [PATCH 11/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/likelihoods/negativebinomial.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index ea47243..66f559e 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -36,7 +36,9 @@ end @functor NegativeBinomialLikelihood function (l::NegativeBinomialLikelihood)(::Real) - error("not implemented for type $(typeof(l)). See `NegativeBinomialLikelihood` docs") + return error( + "not implemented for type $(typeof(l)). See `NegativeBinomialLikelihood` docs" + ) end (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) From b512b69bee9e8e127b9a8b8ba9785f5441df8661 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 19 Apr 2022 10:20:55 +0200 Subject: [PATCH 12/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/likelihoods/negativebinomial.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 66f559e..18fbc0c 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -43,7 +43,6 @@ end (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) - """ NBParamI From 67c4cb39f22755d718b230428ecd9b39f772a2ca Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Tue, 19 Apr 2022 11:14:09 +0200 Subject: [PATCH 13/44] Work around makefunctor --- src/likelihoods/negativebinomial.jl | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 18fbc0c..490b8d4 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -33,14 +33,25 @@ end l; successes ) -@functor NegativeBinomialLikelihood - function (l::NegativeBinomialLikelihood)(::Real) return error( "not implemented for type $(typeof(l)). See `NegativeBinomialLikelihood` docs" ) end +# Workaround for https://github.com/FluxML/Functors.jl/issues/40 +function Functors.makefunctor(m::Module, T::Type{<:Tlik}, fs=fieldnames(T)) where {Tparam,Tlik<:NegativeBinomialLikelihood{Tparam}} + yᵢ = 0 + escargs = map(fieldnames(T)) do f + :(y[$(yᵢ += 1)]) + end + escfs = [:($f = x.$f) for f in fs] + + @eval m begin + $Functors.functor(::Type{<:$T}, x) = ($(escfs...),), y -> $T($(escargs...)) + end +end + (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) """ @@ -102,3 +113,4 @@ end function NegativeBinomialLikelihood{NBParamIII}(l=exp; successes::Real=1) return NegativeBinomialLikelihood{NBParamIII}((; successes), l) end + From baeda7b28ee30697ee679f7b15156f5d95ffff6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 19 Apr 2022 11:18:22 +0200 Subject: [PATCH 14/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/likelihoods/negativebinomial.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 490b8d4..b62dd1b 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -40,7 +40,9 @@ function (l::NegativeBinomialLikelihood)(::Real) end # Workaround for https://github.com/FluxML/Functors.jl/issues/40 -function Functors.makefunctor(m::Module, T::Type{<:Tlik}, fs=fieldnames(T)) where {Tparam,Tlik<:NegativeBinomialLikelihood{Tparam}} +function Functors.makefunctor( + m::Module, T::Type{<:Tlik}, fs=fieldnames(T) +) where {Tparam,Tlik<:NegativeBinomialLikelihood{Tparam}} yᵢ = 0 escargs = map(fieldnames(T)) do f :(y[$(yᵢ += 1)]) From c8160c7bef4daa4353ecbfed82c2c6c5cac07c6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 19 Apr 2022 11:18:27 +0200 Subject: [PATCH 15/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/likelihoods/negativebinomial.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index b62dd1b..38784db 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -115,4 +115,3 @@ end function NegativeBinomialLikelihood{NBParamIII}(l=exp; successes::Real=1) return NegativeBinomialLikelihood{NBParamIII}((; successes), l) end - From 769543380591dcf641b4c07453280a42c2201183 Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Tue, 19 Apr 2022 11:58:12 +0200 Subject: [PATCH 16/44] Fix functor function --- src/likelihoods/negativebinomial.jl | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 38784db..cca6180 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -40,18 +40,13 @@ function (l::NegativeBinomialLikelihood)(::Real) end # Workaround for https://github.com/FluxML/Functors.jl/issues/40 -function Functors.makefunctor( - m::Module, T::Type{<:Tlik}, fs=fieldnames(T) +function Functors.functor( + ::Type{<:Tlik}, x ) where {Tparam,Tlik<:NegativeBinomialLikelihood{Tparam}} - yᵢ = 0 - escargs = map(fieldnames(T)) do f - :(y[$(yᵢ += 1)]) - end - escfs = [:($f = x.$f) for f in fs] - - @eval m begin - $Functors.functor(::Type{<:$T}, x) = ($(escfs...),), y -> $T($(escargs...)) + function reconstruct_lik(xs) + NegativeBinomialLikelihood{Tparam}(xs.params, xs.invlink) end + return (params = x.params, invlink = x.invlink), reconstruct_lik end (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) From 65f37341ea2169af3a1e9ec0c987bf848f83fdb3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 19 Apr 2022 12:13:57 +0200 Subject: [PATCH 17/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/likelihoods/negativebinomial.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index cca6180..8e6a871 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -44,7 +44,7 @@ function Functors.functor( ::Type{<:Tlik}, x ) where {Tparam,Tlik<:NegativeBinomialLikelihood{Tparam}} function reconstruct_lik(xs) - NegativeBinomialLikelihood{Tparam}(xs.params, xs.invlink) + return NegativeBinomialLikelihood{Tparam}(xs.params, xs.invlink) end return (params = x.params, invlink = x.invlink), reconstruct_lik end From 4460c8e134ccbe607f23b4efce3400c8f63ecd0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 19 Apr 2022 12:14:02 +0200 Subject: [PATCH 18/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/likelihoods/negativebinomial.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 8e6a871..2077c20 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -46,7 +46,7 @@ function Functors.functor( function reconstruct_lik(xs) return NegativeBinomialLikelihood{Tparam}(xs.params, xs.invlink) end - return (params = x.params, invlink = x.invlink), reconstruct_lik + return (params=x.params, invlink=x.invlink), reconstruct_lik end (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) From 43766773d1e171a8e0986460c35665f40f8777e2 Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Tue, 19 Apr 2022 12:18:39 +0200 Subject: [PATCH 19/44] Add test for error on custom NBParam --- test/likelihoods/negativebinomial.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index e3da07d..79a6d02 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -3,7 +3,6 @@ ((NBParamI, :successes), (NBParamII, :failures), (NBParamIII, :successes)) @testset "$(nameof(nbparam))" begin @eval begin - sym_kwarg = $(Meta.quot(kwarg)) for args in ((logistic,), (LogisticLink(),)), kwargs in ((), (; $(kwarg)=1)) lik = NegativeBinomialLikelihood{$(nbparam)}(args...; kwargs...) @test lik isa NegativeBinomialLikelihood{ @@ -23,4 +22,9 @@ end end end + struct NBParamFoo <: GPLikelihoods.NBParam end + function NegativeBinomialLikelihood{NBParamFoo}(l=logistic; bar=1) + return NegativeBinomialLikelihood{NBParamFoo}((; bar), GPLikelihoods.link(l)) + end + @test_throws ErrorException NegativeBinomialLikelihood{NBParamFoo}()(2.0) end From 94921fbe09aae38a289592078e795eab657e872f Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Fri, 22 Apr 2022 11:45:36 +0200 Subject: [PATCH 20/44] Refactor the parametrization --- src/likelihoods/negativebinomial.jl | 85 ++++++++++++++-------------- test/likelihoods/negativebinomial.jl | 30 +++------- 2 files changed, 48 insertions(+), 67 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 2077c20..ad890e7 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -1,36 +1,47 @@ abstract type NBParam end """ - NegativeBinomialLikelihood{NBParam}(l=logistic/exp; kwargs...) + NegativeBinomialLikelihood(param::NBParam, l::Function/Link) There are many possible parametrizations for the Negative Binomial likelihood. The `NegativeBinomialLikelihood` has a special structure, the first type `NBParam` -defines what parametrization is used. +defines what parametrization is used, and contains the related parameters. Each `NBParam` has its own documentation: - [`NBParamI`](@ref): This is the definition used by Distributions.jl. - [`NBParamII`](@ref): This is the definition used by Wikipedia. - [`NBParamIII`](@ref): This is the definition based on the mean. To create a new parametrization, you need to; -- Create a new typee from `struct MyNBParam <: NBParam end` -- Dispatch `(l::NegativeBinomialLikelihood{MyNBParam})(f::Real)`, which return the [`NegativeBinomial`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial) from `Distributions.jl`. -`NegativeBinomial` follows the parametrization [`NBParamI`](@ref), i.e. the first argument is the number of successes +- Create a new type `struct MyNBParam{T} <: NBParam; myparams::T; end` +- Dispatch `(l::NegativeBinomialLikelihood{<:MyNBParam})(f::Real)`, which return a [`NegativeBinomial`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial) from `Distributions.jl`. +`NegativeBinomial` follows the parametrization of [`NBParamI`](@ref), i.e. the first argument is the number of successes and the second argument is the probability of success. -- Write a constructor for `NegativeBinomialLikelihood{MyNBParam}(link; kwargs...)` + +## Examples + +```julia-repl +julia> lik = NegativeBinomialLikelihood(NBParamI(10), logistic) +NegativeBinomialLikelihood{NBParamI{Int64}, LogisticLink}(NBParamI{Int64}(10), LogisticLink(LogExpFunctions.logistic)) +julia> lik(2.0) +Distributions.NegativeBinomial{Float64}(r=10.0, p=0.8807970779778824) + +julia> lik = NegativeBinomialLikelihood(NBParamIII(10), exp) +NegativeBinomialLikelihood{NBParamIII{Int64}, ExpLink}(NBParamIII{Int64}(10), ExpLink(exp)) +julia> lik(2.0) +Distributions.NegativeBinomial{Float64}(r=10.0, p=0.4249256576603398) +``` """ -struct NegativeBinomialLikelihood{NBParam,Tl<:AbstractLink,T} <: AbstractLikelihood - params::T # Likelihood parameters (depends of NBParam) +struct NegativeBinomialLikelihood{Tp<:NBParam,Tl<:AbstractLink} <: AbstractLikelihood + params::Tp # Likelihood parametrization (and parameters) invlink::Tl - function NegativeBinomialLikelihood{Tparam}(params, invlink) where {Tparam} + function NegativeBinomialLikelihood(params::Tparam, invlink) where {Tparam} invlink = link(invlink) - return new{Tparam,typeof(invlink),typeof(params)}(params, invlink) + return new{Tparam,typeof(invlink)}(params, invlink) end end -@deprecate NegativeBinomialLikelihood(l=logistic; successes=1) NegativeBinomialLikelihood{ - NBParamI -}( - l; successes +@deprecate NegativeBinomialLikelihood(l=logistic; successes=1) NegativeBinomialLikelihood( + NBParamI(successes), l ) function (l::NegativeBinomialLikelihood)(::Real) @@ -39,20 +50,12 @@ function (l::NegativeBinomialLikelihood)(::Real) ) end -# Workaround for https://github.com/FluxML/Functors.jl/issues/40 -function Functors.functor( - ::Type{<:Tlik}, x -) where {Tparam,Tlik<:NegativeBinomialLikelihood{Tparam}} - function reconstruct_lik(xs) - return NegativeBinomialLikelihood{Tparam}(xs.params, xs.invlink) - end - return (params=x.params, invlink=x.invlink), reconstruct_lik -end +@functor NegativeBinomialLikelihood (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) """ - NBParamI + NBParamI(successes) Negative Binomial parametrization with `successes` the number of successes and `l(f)` the probability of `success`. @@ -62,18 +65,16 @@ This corresponds to the definition used by [Distributions.jl](https://juliastats p(k|successes, f) = \\frac{\\Gamma(k+successes)}{k! \\Gamma(successes)} l(f)^successes (1 - l(f))^k ``` """ -struct NBParamI <: NBParam end - -function (l::NegativeBinomialLikelihood{NBParamI})(f::Real) - return NegativeBinomial(l.params.successes, l.invlink(f)) +struct NBParamI{T} <: NBParam + successes::T end -function NegativeBinomialLikelihood{NBParamI}(l=logistic; successes::Real=1) - return NegativeBinomialLikelihood{NBParamI}((; successes), l) +function (l::NegativeBinomialLikelihood{<:NBParamI})(f::Real) + return NegativeBinomial(l.params.successes, l.invlink(f)) end """ - NBParamII + NBParamII(failures) Negative Binomial parametrization with `failures` the number of failures and `l(f)` the probability of `success`. @@ -83,30 +84,26 @@ This corresponds to the definition used by [Wikipedia](https://en.wikipedia.org/ p(k|failures, f) = \\frac{\\Gamma(k+failures)}{k! \\Gamma(failure)} l(f)^k (1 - l(f))^{failures} ``` """ -struct NBParamII <: NBParam end - -function (l::NegativeBinomialLikelihood{NBParamII})(f::Real) - return NegativeBinomial(l.params.failures, 1 - l.invlink(f)) +struct NBParamII{T} <: NBParam + failures::T end -function NegativeBinomialLikelihood{NBParamII}(l=logistic; failures::Real=1) - return NegativeBinomialLikelihood{NBParamII}((; failures), l) +function (l::NegativeBinomialLikelihood{<:NBParamII})(f::Real) + return NegativeBinomial(l.params.failures, 1 - l.invlink(f)) end """ - NBParamIII + NBParamIII(successes) Negative Binomial parametrization with mean `μ=l(f)` and number of successes `successes`. See the definition given in the [Wikipedia article](https://en.wikipedia.org/wiki/Negative_binomial_distribution#Alternative_formulations) """ -struct NBParamIII <: NBParam end +struct NBParamIII{T} <: NBParam + successes::T +end -function (l::NegativeBinomialLikelihood{NBParamIII})(f::Real) +function (l::NegativeBinomialLikelihood{<:NBParamIII})(f::Real) μ = l.invlink(f) r = l.params.successes return NegativeBinomial(r, μ / (μ + r)) end - -function NegativeBinomialLikelihood{NBParamIII}(l=exp; successes::Real=1) - return NegativeBinomialLikelihood{NBParamIII}((; successes), l) -end diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index 79a6d02..aebd7af 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -1,30 +1,14 @@ @testset "NegativeBinomialLikelihood" begin - for (nbparam, kwarg) in - ((NBParamI, :successes), (NBParamII, :failures), (NBParamIII, :successes)) - @testset "$(nameof(nbparam))" begin - @eval begin - for args in ((logistic,), (LogisticLink(),)), kwargs in ((), (; $(kwarg)=1)) - lik = NegativeBinomialLikelihood{$(nbparam)}(args...; kwargs...) - @test lik isa NegativeBinomialLikelihood{ - $(nbparam),LogisticLink,<:NamedTuple{<:Any,<:Tuple{Int}} - } - end - - for args in ((normcdf,), (NormalCDFLink(),)), kwargs in ((; $(kwarg)=2.0),) - lik = NegativeBinomialLikelihood{$(nbparam)}(args...; kwargs...) - @test lik isa NegativeBinomialLikelihood{ - $(nbparam),NormalCDFLink,<:NamedTuple{<:Any,<:Tuple{Float64}} - } - end - - lik = NegativeBinomialLikelihood{$(nbparam)}() + rs = (10, 9.5) + for nbparam in (NBParamI, NBParamII, NBParamIII) + for r in (10, 9.5) # Test both input types + @testset "$(nameof(nbparam)), r=$r" begin + lik = NegativeBinomialLikelihood(nbparam(r), logistic) + @test lik isa NegativeBinomialLikelihood{<:nbparam} test_interface(lik, NegativeBinomial; functor_args=(:params, :invlink)) end end end struct NBParamFoo <: GPLikelihoods.NBParam end - function NegativeBinomialLikelihood{NBParamFoo}(l=logistic; bar=1) - return NegativeBinomialLikelihood{NBParamFoo}((; bar), GPLikelihoods.link(l)) - end - @test_throws ErrorException NegativeBinomialLikelihood{NBParamFoo}()(2.0) + @test_throws ErrorException NegativeBinomialLikelihood(NBParamFoo(), logistic)(2.0) end From 72136feb99add1236a117db2a3f5bd6a1ec1390d Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Fri, 22 Apr 2022 18:51:28 +0200 Subject: [PATCH 21/44] Readapt the interface again! --- src/GPLikelihoods.jl | 6 +- src/likelihoods/negativebinomial.jl | 101 ++++++++++++++++++++------- test/likelihoods/negativebinomial.jl | 13 +++- 3 files changed, 91 insertions(+), 29 deletions(-) diff --git a/src/GPLikelihoods.jl b/src/GPLikelihoods.jl index 38fce29..ed79c5a 100644 --- a/src/GPLikelihoods.jl +++ b/src/GPLikelihoods.jl @@ -15,7 +15,11 @@ export BernoulliLikelihood, ExponentialLikelihood, GammaLikelihood, NegativeBinomialLikelihood -export NBParamI, NBParamII, NBParamIII +export NBParamI, + NBParamII, + NBParamPower, + NBParamSuccess, + NBParamFailure export Link, ChainLink, BijectiveSimplexLink, diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index ad890e7..bcbd8ec 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -4,12 +4,22 @@ abstract type NBParam end NegativeBinomialLikelihood(param::NBParam, l::Function/Link) There are many possible parametrizations for the Negative Binomial likelihood. +We follow the convention laid out in p.137 of [^1] and some common parametrizations. The `NegativeBinomialLikelihood` has a special structure, the first type `NBParam` defines what parametrization is used, and contains the related parameters. -Each `NBParam` has its own documentation: -- [`NBParamI`](@ref): This is the definition used by Distributions.jl. -- [`NBParamII`](@ref): This is the definition used by Wikipedia. -- [`NBParamIII`](@ref): This is the definition based on the mean. + +## `NBParam` predefined types + +### `NBParam` with `p = link(f)` the probability of success +- [`NBParamSuccess`](@ref): This is the definition used in [`Distributions.jl`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial). +- [`NBParamFailure`](@ref): This is the definition used in [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution). + + +### `NBParam` with `mean = link(f)` +- [`NBParamI`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α)` +- [`NBParamII`](@ref): Mean is linked to `f` and variance is given by `μ(1 + αμ)` +- [`NBParamPower`](@ref): Mean is linked to `f` and variance is given by `μ + αμ^ρ` + To create a new parametrization, you need to; - Create a new type `struct MyNBParam{T} <: NBParam; myparams::T; end` @@ -20,16 +30,20 @@ and the second argument is the probability of success. ## Examples ```julia-repl -julia> lik = NegativeBinomialLikelihood(NBParamI(10), logistic) -NegativeBinomialLikelihood{NBParamI{Int64}, LogisticLink}(NBParamI{Int64}(10), LogisticLink(LogExpFunctions.logistic)) -julia> lik(2.0) -Distributions.NegativeBinomial{Float64}(r=10.0, p=0.8807970779778824) - -julia> lik = NegativeBinomialLikelihood(NBParamIII(10), exp) -NegativeBinomialLikelihood{NBParamIII{Int64}, ExpLink}(NBParamIII{Int64}(10), ExpLink(exp)) -julia> lik(2.0) -Distributions.NegativeBinomial{Float64}(r=10.0, p=0.4249256576603398) +julia> NegativeBinomialLikelihood(NBParamSuccess(10), logistic)(2.0) +NegativeBinomial{Float64}(r=10.0, p=0.8807970779778824) +julia> NegativeBinomialLikelihood(NBParamFailure(10), logistic)(2.0) +NegativeBinomial{Float64}(r=10.0, p=0.11920292202211757) + +julia> d = NegativeBinomialLikelihood(NBParamI(3.0), exp)(2.0) +NegativeBinomial{Float64}(r=2.4630186996435506, p=0.25) +julia> mean(d) ≈ exp(2.0) +true +julia> var(d) ≈ exp(2.0) * (1 + 3.0) +true ``` + +[^1] Hilbe, Joseph M. Negative binomial regression. Cambridge University Press, 2011. """ struct NegativeBinomialLikelihood{Tp<:NBParam,Tl<:AbstractLink} <: AbstractLikelihood params::Tp # Likelihood parametrization (and parameters) @@ -55,7 +69,7 @@ end (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) """ - NBParamI(successes) + NBParamSuccess(successes) Negative Binomial parametrization with `successes` the number of successes and `l(f)` the probability of `success`. @@ -65,16 +79,16 @@ This corresponds to the definition used by [Distributions.jl](https://juliastats p(k|successes, f) = \\frac{\\Gamma(k+successes)}{k! \\Gamma(successes)} l(f)^successes (1 - l(f))^k ``` """ -struct NBParamI{T} <: NBParam +struct NBParamSuccess{T} <: NBParam successes::T end -function (l::NegativeBinomialLikelihood{<:NBParamI})(f::Real) +function (l::NegativeBinomialLikelihood{<:NBParamSuccess})(f::Real) return NegativeBinomial(l.params.successes, l.invlink(f)) end """ - NBParamII(failures) + NBParamFailure(failures) Negative Binomial parametrization with `failures` the number of failures and `l(f)` the probability of `success`. @@ -84,26 +98,59 @@ This corresponds to the definition used by [Wikipedia](https://en.wikipedia.org/ p(k|failures, f) = \\frac{\\Gamma(k+failures)}{k! \\Gamma(failure)} l(f)^k (1 - l(f))^{failures} ``` """ -struct NBParamII{T} <: NBParam +struct NBParamFailure{T} <: NBParam failures::T end -function (l::NegativeBinomialLikelihood{<:NBParamII})(f::Real) +function (l::NegativeBinomialLikelihood{<:NBParamFailure})(f::Real) return NegativeBinomial(l.params.failures, 1 - l.invlink(f)) end +# Helper function to convert mean and variance to p and r +_nb_mean_var_to_r_p(μ::Real, v::Real) = abs2(μ) / (v - μ), μ / v + """ - NBParamIII(successes) + NBParamI(α) -Negative Binomial parametrization with mean `μ=l(f)` and number of successes `successes`. -See the definition given in the [Wikipedia article](https://en.wikipedia.org/wiki/Negative_binomial_distribution#Alternative_formulations) +Negative Binomial parametrization with mean `μ=l(f)` and variance `v=μ(1 + α)`. """ -struct NBParamIII{T} <: NBParam - successes::T +struct NBParamI{T} <: NBParam + α::T +end + +function (l::NegativeBinomialLikelihood{<:NBParamI})(f::Real) + μ = l.invlink(f) + v = μ * (1 + l.params.α) + return NegativeBinomial(_nb_mean_var_to_r_p(μ, v)...) +end + +""" + NBParamII(α) + +Negative Binomial parametrization with mean `μ=l(f)` and variance `v=μ(1 + αμ)`. +""" +struct NBParamII{T} <: NBParam + α::T +end + +function (l::NegativeBinomialLikelihood{<:NBParamII})(f::Real) + μ = l.invlink(f) + v = μ * (1 + l.params.α * μ) + return NegativeBinomial(_nb_mean_var_to_r_p(μ, v)...) +end + +""" + NBParamPower(α, ρ) + +Negative Binomial parametrization with mean `μ=l(f)` and variance `v=μ(1 + αμ^ρ)`. +""" +struct NBParamPower{Tα,Tρ} <: NBParam + α::Tα + ρ::Tρ end -function (l::NegativeBinomialLikelihood{<:NBParamIII})(f::Real) +function (l::NegativeBinomialLikelihood{<:NBParamPower})(f::Real) μ = l.invlink(f) - r = l.params.successes - return NegativeBinomial(r, μ / (μ + r)) + v = μ * (1 + l.params.α * μ^l.params.ρ) + return NegativeBinomial(_nb_mean_var_to_r_p(μ, v)...) end diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index aebd7af..b50aae4 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -1,6 +1,7 @@ @testset "NegativeBinomialLikelihood" begin rs = (10, 9.5) - for nbparam in (NBParamI, NBParamII, NBParamIII) + # Test based on p success + for nbparam in (NBParamSuccess, NBParamFailure) for r in (10, 9.5) # Test both input types @testset "$(nameof(nbparam)), r=$r" begin lik = NegativeBinomialLikelihood(nbparam(r), logistic) @@ -9,6 +10,16 @@ end end end + # Test based on mean = link(f) + for (nbparam, args) in ((NBParamI, (2.0,)), (NBParamII, (3.0,)), (NBParamPower, (2.0, 2.0))) + @testset "$(nameof(nbparam))" begin + lik = NegativeBinomialLikelihood(nbparam(args...), exp) + @test lik isa NegativeBinomialLikelihood{<:nbparam} + x = rand() + @test mean(lik(x)) ≈ exp(x) + test_interface(lik, NegativeBinomial; functor_args=(:params, :invlink)) + end + end struct NBParamFoo <: GPLikelihoods.NBParam end @test_throws ErrorException NegativeBinomialLikelihood(NBParamFoo(), logistic)(2.0) end From 654c384cb862bed043e92c860257064fed84cf4f Mon Sep 17 00:00:00 2001 From: Theo Galy-Fajou Date: Fri, 22 Apr 2022 18:53:54 +0200 Subject: [PATCH 22/44] Formatting --- src/GPLikelihoods.jl | 6 +----- test/likelihoods/negativebinomial.jl | 3 ++- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/src/GPLikelihoods.jl b/src/GPLikelihoods.jl index ed79c5a..c385e71 100644 --- a/src/GPLikelihoods.jl +++ b/src/GPLikelihoods.jl @@ -15,11 +15,7 @@ export BernoulliLikelihood, ExponentialLikelihood, GammaLikelihood, NegativeBinomialLikelihood -export NBParamI, - NBParamII, - NBParamPower, - NBParamSuccess, - NBParamFailure +export NBParamI, NBParamII, NBParamPower, NBParamSuccess, NBParamFailure export Link, ChainLink, BijectiveSimplexLink, diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index b50aae4..6ab5057 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -11,7 +11,8 @@ end end # Test based on mean = link(f) - for (nbparam, args) in ((NBParamI, (2.0,)), (NBParamII, (3.0,)), (NBParamPower, (2.0, 2.0))) + for (nbparam, args) in + ((NBParamI, (2.0,)), (NBParamII, (3.0,)), (NBParamPower, (2.0, 2.0))) @testset "$(nameof(nbparam))" begin lik = NegativeBinomialLikelihood(nbparam(args...), exp) @test lik isa NegativeBinomialLikelihood{<:nbparam} From 2b0385c9ac6d225f657ac776c7b850ba1c6aaa2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Mon, 25 Apr 2022 14:23:50 +0200 Subject: [PATCH 23/44] Add additional subtypes --- src/likelihoods/negativebinomial.jl | 25 +++++++++++++++++-------- test/likelihoods/negativebinomial.jl | 11 +++++++---- 2 files changed, 24 insertions(+), 12 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index bcbd8ec..d0ec9b6 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -1,4 +1,10 @@ abstract type NBParam end +default_invlink(::NBParam) = error("Not implemented") + +abstract type NBParamProb <: NBParam end +default_invlink(::NBParamProb) = logistic +abstract type NBParamMean <: NBParam end +default_invlink(::NBParamMean) = exp """ NegativeBinomialLikelihood(param::NBParam, l::Function/Link) @@ -7,15 +13,18 @@ There are many possible parametrizations for the Negative Binomial likelihood. We follow the convention laid out in p.137 of [^1] and some common parametrizations. The `NegativeBinomialLikelihood` has a special structure, the first type `NBParam` defines what parametrization is used, and contains the related parameters. +`NBParam` itself has two subtypes: +- `NBParamProb` for parametrizations where `f->p`, the probability of success +- `NBParamMean` for parametrizations where `f->μ`, the mean ## `NBParam` predefined types -### `NBParam` with `p = link(f)` the probability of success +### `NBParamProb` types with `p = link(f)` the probability of success - [`NBParamSuccess`](@ref): This is the definition used in [`Distributions.jl`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial). - [`NBParamFailure`](@ref): This is the definition used in [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution). -### `NBParam` with `mean = link(f)` +### `NBParamMean` types with `mean = link(f)` - [`NBParamI`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α)` - [`NBParamII`](@ref): Mean is linked to `f` and variance is given by `μ(1 + αμ)` - [`NBParamPower`](@ref): Mean is linked to `f` and variance is given by `μ + αμ^ρ` @@ -48,7 +57,7 @@ true struct NegativeBinomialLikelihood{Tp<:NBParam,Tl<:AbstractLink} <: AbstractLikelihood params::Tp # Likelihood parametrization (and parameters) invlink::Tl - function NegativeBinomialLikelihood(params::Tparam, invlink) where {Tparam} + function NegativeBinomialLikelihood(params::Tparam, invlink=default_invlink(params)) where {Tparam<:NBParam} invlink = link(invlink) return new{Tparam,typeof(invlink)}(params, invlink) end @@ -79,7 +88,7 @@ This corresponds to the definition used by [Distributions.jl](https://juliastats p(k|successes, f) = \\frac{\\Gamma(k+successes)}{k! \\Gamma(successes)} l(f)^successes (1 - l(f))^k ``` """ -struct NBParamSuccess{T} <: NBParam +struct NBParamSuccess{T} <: NBParamProb successes::T end @@ -98,7 +107,7 @@ This corresponds to the definition used by [Wikipedia](https://en.wikipedia.org/ p(k|failures, f) = \\frac{\\Gamma(k+failures)}{k! \\Gamma(failure)} l(f)^k (1 - l(f))^{failures} ``` """ -struct NBParamFailure{T} <: NBParam +struct NBParamFailure{T} <: NBParamProb failures::T end @@ -114,7 +123,7 @@ _nb_mean_var_to_r_p(μ::Real, v::Real) = abs2(μ) / (v - μ), μ / v Negative Binomial parametrization with mean `μ=l(f)` and variance `v=μ(1 + α)`. """ -struct NBParamI{T} <: NBParam +struct NBParamI{T} <: NBParamMean α::T end @@ -129,7 +138,7 @@ end Negative Binomial parametrization with mean `μ=l(f)` and variance `v=μ(1 + αμ)`. """ -struct NBParamII{T} <: NBParam +struct NBParamII{T} <: NBParamMean α::T end @@ -144,7 +153,7 @@ end Negative Binomial parametrization with mean `μ=l(f)` and variance `v=μ(1 + αμ^ρ)`. """ -struct NBParamPower{Tα,Tρ} <: NBParam +struct NBParamPower{Tα,Tρ} <: NBParamMean α::Tα ρ::Tρ end diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index 6ab5057..5833eb0 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -5,17 +5,20 @@ for r in (10, 9.5) # Test both input types @testset "$(nameof(nbparam)), r=$r" begin lik = NegativeBinomialLikelihood(nbparam(r), logistic) - @test lik isa NegativeBinomialLikelihood{<:nbparam} + @test lik isa NegativeBinomialLikelihood{<:nbparam,LogisticLink} + lik = NegativeBinomialLikelihood(nbparam(args...)) + @test lik isa NegativeBinomialLikelihood{<:nbparam,LogisticLink} test_interface(lik, NegativeBinomial; functor_args=(:params, :invlink)) end end end # Test based on mean = link(f) - for (nbparam, args) in - ((NBParamI, (2.0,)), (NBParamII, (3.0,)), (NBParamPower, (2.0, 2.0))) + for (nbparam, args) in ((NBParamI, (2.0,)), (NBParamII, (3.0,)), (NBParamPower, (2.0, 2.0))) @testset "$(nameof(nbparam))" begin lik = NegativeBinomialLikelihood(nbparam(args...), exp) - @test lik isa NegativeBinomialLikelihood{<:nbparam} + @test lik isa NegativeBinomialLikelihood{<:nbparam,ExpLink} + lik = NegativeBinomialLikelihood(nbparam(args...)) + @test lik isa NegativeBinomialLikelihood{<:nbparam,ExpLink} x = rand() @test mean(lik(x)) ≈ exp(x) test_interface(lik, NegativeBinomial; functor_args=(:params, :invlink)) From 4ed2a5d4a83c681e2b34c85f58daf346db8dfbf2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Mon, 25 Apr 2022 14:26:52 +0200 Subject: [PATCH 24/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/likelihoods/negativebinomial.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index d0ec9b6..81647dc 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -57,7 +57,9 @@ true struct NegativeBinomialLikelihood{Tp<:NBParam,Tl<:AbstractLink} <: AbstractLikelihood params::Tp # Likelihood parametrization (and parameters) invlink::Tl - function NegativeBinomialLikelihood(params::Tparam, invlink=default_invlink(params)) where {Tparam<:NBParam} + function NegativeBinomialLikelihood( + params::Tparam, invlink=default_invlink(params) + ) where {Tparam<:NBParam} invlink = link(invlink) return new{Tparam,typeof(invlink)}(params, invlink) end From c6bf273281808a0ee07ed03a8e661eba78ae788c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Mon, 25 Apr 2022 14:32:34 +0200 Subject: [PATCH 25/44] Fix and format --- src/likelihoods/negativebinomial.jl | 4 +++- test/likelihoods/negativebinomial.jl | 11 ++++++----- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index d0ec9b6..81647dc 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -57,7 +57,9 @@ true struct NegativeBinomialLikelihood{Tp<:NBParam,Tl<:AbstractLink} <: AbstractLikelihood params::Tp # Likelihood parametrization (and parameters) invlink::Tl - function NegativeBinomialLikelihood(params::Tparam, invlink=default_invlink(params)) where {Tparam<:NBParam} + function NegativeBinomialLikelihood( + params::Tparam, invlink=default_invlink(params) + ) where {Tparam<:NBParam} invlink = link(invlink) return new{Tparam,typeof(invlink)}(params, invlink) end diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index 5833eb0..a3f2e87 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -4,21 +4,22 @@ for nbparam in (NBParamSuccess, NBParamFailure) for r in (10, 9.5) # Test both input types @testset "$(nameof(nbparam)), r=$r" begin - lik = NegativeBinomialLikelihood(nbparam(r), logistic) + lik = NegativeBinomialLikelihood(nbparam(r)) @test lik isa NegativeBinomialLikelihood{<:nbparam,LogisticLink} - lik = NegativeBinomialLikelihood(nbparam(args...)) + lik = NegativeBinomialLikelihood(nbparam(r), logistic) @test lik isa NegativeBinomialLikelihood{<:nbparam,LogisticLink} test_interface(lik, NegativeBinomial; functor_args=(:params, :invlink)) end end end # Test based on mean = link(f) - for (nbparam, args) in ((NBParamI, (2.0,)), (NBParamII, (3.0,)), (NBParamPower, (2.0, 2.0))) + for (nbparam, args) in + ((NBParamI, (2.0,)), (NBParamII, (3.0,)), (NBParamPower, (2.0, 2.0))) @testset "$(nameof(nbparam))" begin - lik = NegativeBinomialLikelihood(nbparam(args...), exp) - @test lik isa NegativeBinomialLikelihood{<:nbparam,ExpLink} lik = NegativeBinomialLikelihood(nbparam(args...)) @test lik isa NegativeBinomialLikelihood{<:nbparam,ExpLink} + lik = NegativeBinomialLikelihood(nbparam(args...), exp) + @test lik isa NegativeBinomialLikelihood{<:nbparam,ExpLink} x = rand() @test mean(lik(x)) ≈ exp(x) test_interface(lik, NegativeBinomial; functor_args=(:params, :invlink)) From e0b48de3cfe6232f7fc5f18a1c1a639e68bb0d61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Mon, 25 Apr 2022 15:34:05 +0200 Subject: [PATCH 26/44] Adds more docs --- src/likelihoods/negativebinomial.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 81647dc..5f88359 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -60,6 +60,7 @@ struct NegativeBinomialLikelihood{Tp<:NBParam,Tl<:AbstractLink} <: AbstractLikel function NegativeBinomialLikelihood( params::Tparam, invlink=default_invlink(params) ) where {Tparam<:NBParam} + # we convert `invlink` into `Link` if it's not the case already invlink = link(invlink) return new{Tparam,typeof(invlink)}(params, invlink) end @@ -71,7 +72,9 @@ end function (l::NegativeBinomialLikelihood)(::Real) return error( - "not implemented for type $(typeof(l)). See `NegativeBinomialLikelihood` docs" + "not implemented for type $(typeof(l)). For your custom type to run you ", + "need to implement `(l::NegativeBinomialLikelihood{<:MyNBParam})(f::Real)`.", + "For a full explanation, see `NegativeBinomialLikelihood` docs", ) end From 393198b8cbe51a5c51cffe885e79193ef27ca362 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 26 Apr 2022 21:06:10 +0200 Subject: [PATCH 27/44] Update api.md --- docs/src/api.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/src/api.md b/docs/src/api.md index 1890a4d..e48ccd5 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -19,9 +19,11 @@ PoissonLikelihood ```@docs NegativeBinomialLikelihood +NBParamSuccess +NBParamFailure NBParamI NBParamII -NBParamIII +NBParamPower ``` ## Links @@ -51,4 +53,4 @@ LogisticLink ProbitLink NormalCDFLink SoftMaxLink -``` \ No newline at end of file +``` From 49e2e1760b53e3f0ad4ce7e587b8d0aa0141ad47 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 26 Apr 2022 21:07:54 +0200 Subject: [PATCH 28/44] Apply suggestions from code review Co-authored-by: st-- --- src/likelihoods/negativebinomial.jl | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 5f88359..645a060 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -1,5 +1,4 @@ abstract type NBParam end -default_invlink(::NBParam) = error("Not implemented") abstract type NBParamProb <: NBParam end default_invlink(::NBParamProb) = logistic @@ -7,15 +6,15 @@ abstract type NBParamMean <: NBParam end default_invlink(::NBParamMean) = exp """ - NegativeBinomialLikelihood(param::NBParam, l::Function/Link) + NegativeBinomialLikelihood(param::NBParam, link::Union{Function,Link}) There are many possible parametrizations for the Negative Binomial likelihood. -We follow the convention laid out in p.137 of [^1] and some common parametrizations. -The `NegativeBinomialLikelihood` has a special structure, the first type `NBParam` -defines what parametrization is used, and contains the related parameters. +We follow the convention laid out in p.137 of [^1] and provide some common parametrizations. +The `NegativeBinomialLikelihood` has a special structure; the first type parameter `NBParam` +defines in what parametrization the latent function is used, and contains the other (scalar) parameters. `NBParam` itself has two subtypes: -- `NBParamProb` for parametrizations where `f->p`, the probability of success -- `NBParamMean` for parametrizations where `f->μ`, the mean +- `NBParamProb` for parametrizations where `f -> p`, the probability of success of a Bernoulli event +- `NBParamMean` for parametrizations where `f -> μ`, the mean of the number of events ## `NBParam` predefined types @@ -24,16 +23,16 @@ defines what parametrization is used, and contains the related parameters. - [`NBParamFailure`](@ref): This is the definition used in [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution). -### `NBParamMean` types with `mean = link(f)` +### `NBParamMean` types with `μ = link(f)` the mean/expected number of events - [`NBParamI`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α)` - [`NBParamII`](@ref): Mean is linked to `f` and variance is given by `μ(1 + αμ)` -- [`NBParamPower`](@ref): Mean is linked to `f` and variance is given by `μ + αμ^ρ` +- [`NBParamPower`](@ref): Mean is linked to `f` and variance is given by `μ + α*μ^ρ` -To create a new parametrization, you need to; -- Create a new type `struct MyNBParam{T} <: NBParam; myparams::T; end` -- Dispatch `(l::NegativeBinomialLikelihood{<:MyNBParam})(f::Real)`, which return a [`NegativeBinomial`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial) from `Distributions.jl`. -`NegativeBinomial` follows the parametrization of [`NBParamI`](@ref), i.e. the first argument is the number of successes +To create a new parametrization, you need to: +- create a new type `struct MyNBParam{T} <: NBParam; myparams::T; end`; +- dispatch `(l::NegativeBinomialLikelihood{<:MyNBParam})(f::Real)`, which must return a [`NegativeBinomial`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial) from `Distributions.jl`. +`NegativeBinomial` follows the parametrization of [`NBParamSuccess`](@ref), i.e. the first argument is the number of successes and the second argument is the probability of success. ## Examples @@ -73,7 +72,7 @@ end function (l::NegativeBinomialLikelihood)(::Real) return error( "not implemented for type $(typeof(l)). For your custom type to run you ", - "need to implement `(l::NegativeBinomialLikelihood{<:MyNBParam})(f::Real)`.", + "need to implement `(l::NegativeBinomialLikelihood{<:MyNBParam})(f::Real)`. ", "For a full explanation, see `NegativeBinomialLikelihood` docs", ) end @@ -86,7 +85,7 @@ end NBParamSuccess(successes) Negative Binomial parametrization with `successes` the number of successes and -`l(f)` the probability of `success`. +`link(f)` the probability of `success`. This corresponds to the definition used by [Distributions.jl](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial). ```math From 38009bf65ea8502d9dce1c6b360f28b9cf8064d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 26 Apr 2022 21:23:19 +0200 Subject: [PATCH 29/44] Apply suggestions from code review Co-authored-by: st-- --- src/likelihoods/negativebinomial.jl | 24 ++++++++++++------------ test/likelihoods/negativebinomial.jl | 3 +-- 2 files changed, 13 insertions(+), 14 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 645a060..2b15aec 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -85,11 +85,11 @@ end NBParamSuccess(successes) Negative Binomial parametrization with `successes` the number of successes and -`link(f)` the probability of `success`. +`link(f)` the probability of success. This corresponds to the definition used by [Distributions.jl](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial). ```math - p(k|successes, f) = \\frac{\\Gamma(k+successes)}{k! \\Gamma(successes)} l(f)^successes (1 - l(f))^k + p(k|\\text{successes}, f) = \\frac{\\Gamma(k+\\text{successes})}{k! \\Gamma(\\text{successes})} \\text{link}(f)^\\text{successes} (1 - \text{link}(f))^k ``` """ struct NBParamSuccess{T} <: NBParamProb @@ -104,11 +104,11 @@ end NBParamFailure(failures) Negative Binomial parametrization with `failures` the number of failures and -`l(f)` the probability of `success`. +`link(f)` the probability of success. This corresponds to the definition used by [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution). ```math - p(k|failures, f) = \\frac{\\Gamma(k+failures)}{k! \\Gamma(failure)} l(f)^k (1 - l(f))^{failures} + p(k|failures, f) = \\frac{\\Gamma(k+failures)}{k! \\Gamma(failures)} l(f)^k (1 - l(f))^{failures} ``` """ struct NBParamFailure{T} <: NBParamProb @@ -120,12 +120,12 @@ function (l::NegativeBinomialLikelihood{<:NBParamFailure})(f::Real) end # Helper function to convert mean and variance to p and r -_nb_mean_var_to_r_p(μ::Real, v::Real) = abs2(μ) / (v - μ), μ / v +_nb_mean_excessvar_to_r_p(μ::Real, ev::Real) = μ / ev, 1 / (1 + ev) """ NBParamI(α) -Negative Binomial parametrization with mean `μ=l(f)` and variance `v=μ(1 + α)`. +Negative Binomial parametrization with mean `μ=invlink(f)` and variance `v=μ(1 + α)`. """ struct NBParamI{T} <: NBParamMean α::T @@ -133,8 +133,8 @@ end function (l::NegativeBinomialLikelihood{<:NBParamI})(f::Real) μ = l.invlink(f) - v = μ * (1 + l.params.α) - return NegativeBinomial(_nb_mean_var_to_r_p(μ, v)...) + ev = l.params.α + return NegativeBinomial(_nb_mean_excessvar_to_r_p(μ, ev)...) end """ @@ -148,8 +148,8 @@ end function (l::NegativeBinomialLikelihood{<:NBParamII})(f::Real) μ = l.invlink(f) - v = μ * (1 + l.params.α * μ) - return NegativeBinomial(_nb_mean_var_to_r_p(μ, v)...) + ev = l.params.α * μ + return NegativeBinomial(_nb_mean_excessvar_to_r_p(μ, ev)...) end """ @@ -164,6 +164,6 @@ end function (l::NegativeBinomialLikelihood{<:NBParamPower})(f::Real) μ = l.invlink(f) - v = μ * (1 + l.params.α * μ^l.params.ρ) - return NegativeBinomial(_nb_mean_var_to_r_p(μ, v)...) + ev = l.params.α * μ^l.params.ρ + return NegativeBinomial(_nb_mean_excessvar_to_r_p(μ, ev)...) end diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index a3f2e87..1e55a4d 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -1,5 +1,4 @@ @testset "NegativeBinomialLikelihood" begin - rs = (10, 9.5) # Test based on p success for nbparam in (NBParamSuccess, NBParamFailure) for r in (10, 9.5) # Test both input types @@ -14,7 +13,7 @@ end # Test based on mean = link(f) for (nbparam, args) in - ((NBParamI, (2.0,)), (NBParamII, (3.0,)), (NBParamPower, (2.0, 2.0))) + ((NBParamI, (2.345,)), (NBParamII, (3.456,)), (NBParamPower, (2.34, 3.21))) @testset "$(nameof(nbparam))" begin lik = NegativeBinomialLikelihood(nbparam(args...)) @test lik isa NegativeBinomialLikelihood{<:nbparam,ExpLink} From 79a02d2b54e2478422e40ffaee28c76d89b2a538 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 26 Apr 2022 21:38:37 +0200 Subject: [PATCH 30/44] Update negativebinomial.jl --- src/likelihoods/negativebinomial.jl | 30 ++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 2b15aec..a34fb0e 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -6,7 +6,7 @@ abstract type NBParamMean <: NBParam end default_invlink(::NBParamMean) = exp """ - NegativeBinomialLikelihood(param::NBParam, link::Union{Function,Link}) + NegativeBinomialLikelihood(param::NBParam, invlink::Union{Function,Link}) There are many possible parametrizations for the Negative Binomial likelihood. We follow the convention laid out in p.137 of [^1] and provide some common parametrizations. @@ -14,19 +14,19 @@ The `NegativeBinomialLikelihood` has a special structure; the first type paramet defines in what parametrization the latent function is used, and contains the other (scalar) parameters. `NBParam` itself has two subtypes: - `NBParamProb` for parametrizations where `f -> p`, the probability of success of a Bernoulli event -- `NBParamMean` for parametrizations where `f -> μ`, the mean of the number of events +- `NBParamMean` for parametrizations where `f -> μ`, the expected number of events ## `NBParam` predefined types -### `NBParamProb` types with `p = link(f)` the probability of success +### `NBParamProb` types with `p = invlink(f)` the probability of success - [`NBParamSuccess`](@ref): This is the definition used in [`Distributions.jl`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial). - [`NBParamFailure`](@ref): This is the definition used in [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution). -### `NBParamMean` types with `μ = link(f)` the mean/expected number of events +### `NBParamMean` types with `μ = invlink(f)` the mean/expected number of events - [`NBParamI`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α)` -- [`NBParamII`](@ref): Mean is linked to `f` and variance is given by `μ(1 + αμ)` -- [`NBParamPower`](@ref): Mean is linked to `f` and variance is given by `μ + α*μ^ρ` +- [`NBParamII`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α * μ)` +- [`NBParamPower`](@ref): Mean is linked to `f` and variance is given by `μ + α * μ^ρ` To create a new parametrization, you need to: @@ -59,7 +59,7 @@ struct NegativeBinomialLikelihood{Tp<:NBParam,Tl<:AbstractLink} <: AbstractLikel function NegativeBinomialLikelihood( params::Tparam, invlink=default_invlink(params) ) where {Tparam<:NBParam} - # we convert `invlink` into `Link` if it's not the case already + # we convert `invlink` into a `Link` object if it's not the case already invlink = link(invlink) return new{Tparam,typeof(invlink)}(params, invlink) end @@ -81,15 +81,15 @@ end (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs)) -""" +@doc raw""" NBParamSuccess(successes) Negative Binomial parametrization with `successes` the number of successes and -`link(f)` the probability of success. +`invlink(f)` the probability of success. This corresponds to the definition used by [Distributions.jl](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial). ```math - p(k|\\text{successes}, f) = \\frac{\\Gamma(k+\\text{successes})}{k! \\Gamma(\\text{successes})} \\text{link}(f)^\\text{successes} (1 - \text{link}(f))^k + p(y|\text{successes}, p=invlink(f)) = \frac{\Gamma(y+\text{successes})}{y! \Gamma(\text{successes})} p^\text{successes} (1 - p)^y ``` """ struct NBParamSuccess{T} <: NBParamProb @@ -100,15 +100,15 @@ function (l::NegativeBinomialLikelihood{<:NBParamSuccess})(f::Real) return NegativeBinomial(l.params.successes, l.invlink(f)) end -""" +@doc raw""" NBParamFailure(failures) Negative Binomial parametrization with `failures` the number of failures and -`link(f)` the probability of success. +`invlink(f)` the probability of success. This corresponds to the definition used by [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution). ```math - p(k|failures, f) = \\frac{\\Gamma(k+failures)}{k! \\Gamma(failures)} l(f)^k (1 - l(f))^{failures} + p(y|\text{failures}, p=\text{invlink}(f)) = \frac{\Gamma(y+\text{failures})}{y! \Gamma(\text{failures})} p^y (1 - p)^{\text{failures}} ``` """ struct NBParamFailure{T} <: NBParamProb @@ -140,7 +140,7 @@ end """ NBParamII(α) -Negative Binomial parametrization with mean `μ=l(f)` and variance `v=μ(1 + αμ)`. +Negative Binomial parametrization with mean `μ=invlink(f)` and variance `v=μ(1 + α * μ)`. """ struct NBParamII{T} <: NBParamMean α::T @@ -155,7 +155,7 @@ end """ NBParamPower(α, ρ) -Negative Binomial parametrization with mean `μ=l(f)` and variance `v=μ(1 + αμ^ρ)`. +Negative Binomial parametrization with mean `μ=invlink(f)` and variance `v=μ(1 + α * μ^ρ)`. """ struct NBParamPower{Tα,Tρ} <: NBParamMean α::Tα From 58ac092f8b75545309f147981e1fb4f1fa884731 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 26 Apr 2022 22:32:39 +0200 Subject: [PATCH 31/44] Update test/likelihoods/negativebinomial.jl Co-authored-by: st-- --- test/likelihoods/negativebinomial.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/likelihoods/negativebinomial.jl b/test/likelihoods/negativebinomial.jl index 1e55a4d..ad7f95d 100644 --- a/test/likelihoods/negativebinomial.jl +++ b/test/likelihoods/negativebinomial.jl @@ -19,7 +19,7 @@ @test lik isa NegativeBinomialLikelihood{<:nbparam,ExpLink} lik = NegativeBinomialLikelihood(nbparam(args...), exp) @test lik isa NegativeBinomialLikelihood{<:nbparam,ExpLink} - x = rand() + x = randn() @test mean(lik(x)) ≈ exp(x) test_interface(lik, NegativeBinomial; functor_args=(:params, :invlink)) end From 91212727711ef9676f9f5e8bfada1bc4d6bac105 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 26 Apr 2022 22:52:23 +0200 Subject: [PATCH 32/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: st-- --- src/likelihoods/negativebinomial.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index a34fb0e..612dbac 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -26,7 +26,7 @@ defines in what parametrization the latent function is used, and contains the ot ### `NBParamMean` types with `μ = invlink(f)` the mean/expected number of events - [`NBParamI`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α)` - [`NBParamII`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α * μ)` -- [`NBParamPower`](@ref): Mean is linked to `f` and variance is given by `μ + α * μ^ρ` +- [`NBParamPower`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α * μ^ρ)` To create a new parametrization, you need to: From 739bd3f897785620075526ccc59eab578bab2f69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 26 Apr 2022 22:53:06 +0200 Subject: [PATCH 33/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: st-- --- src/likelihoods/negativebinomial.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 612dbac..c7495da 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -155,7 +155,7 @@ end """ NBParamPower(α, ρ) -Negative Binomial parametrization with mean `μ=invlink(f)` and variance `v=μ(1 + α * μ^ρ)`. +Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α * μ^ρ)`. """ struct NBParamPower{Tα,Tρ} <: NBParamMean α::Tα From db97ebef5f37aca5ab8c8f85993864a4974a5e8d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 26 Apr 2022 22:53:28 +0200 Subject: [PATCH 34/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: st-- --- src/likelihoods/negativebinomial.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index c7495da..75efdd1 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -89,7 +89,7 @@ Negative Binomial parametrization with `successes` the number of successes and This corresponds to the definition used by [Distributions.jl](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial). ```math - p(y|\text{successes}, p=invlink(f)) = \frac{\Gamma(y+\text{successes})}{y! \Gamma(\text{successes})} p^\text{successes} (1 - p)^y + p(y|\text{successes}, p=\text{invlink}(f)) = \frac{\Gamma(y+\text{successes})}{y! \Gamma(\text{successes})} p^\text{successes} (1 - p)^y ``` """ struct NBParamSuccess{T} <: NBParamProb From 5b57eb838571cf9afe6e5c187dc249ffb60527af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 26 Apr 2022 22:53:46 +0200 Subject: [PATCH 35/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: st-- --- src/likelihoods/negativebinomial.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 75efdd1..ea0e569 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -125,7 +125,7 @@ _nb_mean_excessvar_to_r_p(μ::Real, ev::Real) = μ / ev, 1 / (1 + ev) """ NBParamI(α) -Negative Binomial parametrization with mean `μ=invlink(f)` and variance `v=μ(1 + α)`. +Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α)`. """ struct NBParamI{T} <: NBParamMean α::T From f38bd5aa23e7ddf4551b37ec82fe3b609a0855e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 26 Apr 2022 22:54:01 +0200 Subject: [PATCH 36/44] Update src/likelihoods/negativebinomial.jl Co-authored-by: st-- --- src/likelihoods/negativebinomial.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index ea0e569..cfbd9f9 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -140,7 +140,7 @@ end """ NBParamII(α) -Negative Binomial parametrization with mean `μ=invlink(f)` and variance `v=μ(1 + α * μ)`. +Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α * μ)`. """ struct NBParamII{T} <: NBParamMean α::T From 022c891e68aa43351dd74f8055210ee6b76a7094 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Thu, 28 Apr 2022 10:18:59 +0200 Subject: [PATCH 37/44] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f5b9c6d..408e4ce 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GPLikelihoods" uuid = "6031954c-0455-49d7-b3b9-3e1c99afaf40" authors = ["JuliaGaussianProcesses Team"] -version = "0.4.1" +version = "0.4.2" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From 4c8b4d7a3325cc763a0af471d109a5cddc0c4ed3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Thu, 28 Apr 2022 10:21:06 +0200 Subject: [PATCH 38/44] Update negativebinomial.jl --- src/likelihoods/negativebinomial.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index cfbd9f9..685528e 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -125,7 +125,7 @@ _nb_mean_excessvar_to_r_p(μ::Real, ev::Real) = μ / ev, 1 / (1 + ev) """ NBParamI(α) -Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α)`. +Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α)` where `α > 0`. """ struct NBParamI{T} <: NBParamMean α::T @@ -140,7 +140,7 @@ end """ NBParamII(α) -Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α * μ)`. +Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α * μ)` where `α > 0`. """ struct NBParamII{T} <: NBParamMean α::T @@ -155,7 +155,7 @@ end """ NBParamPower(α, ρ) -Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α * μ^ρ)`. +Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α * μ^ρ)` where `α > 0` and `ρ > 0`. """ struct NBParamPower{Tα,Tρ} <: NBParamMean α::Tα From ae77d2868616ec934db40b7d92ec99c29bdb91b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Fri, 29 Apr 2022 08:35:47 +0200 Subject: [PATCH 39/44] Fix different docs issues (#60) * Fix different docs issues * Removed unexisting AbstractLikelihood docs Co-authored-by: st-- --- docs/Project.toml | 6 ++++++ docs/src/api.md | 4 ++-- docs/src/index.md | 28 +++++++++++++++++----------- src/links.jl | 5 +++++ 4 files changed, 30 insertions(+), 13 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 15a2f62..bdd7e95 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,11 @@ [deps] +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +GPLikelihoods = "6031954c-0455-49d7-b3b9-3e1c99afaf40" +StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" [compat] +Distributions = "0.25" Documenter = "0.25, 0.26, 0.27" +GPLikelihoods = "0.4" +StatsFuns = "0.9" diff --git a/docs/src/api.md b/docs/src/api.md index e48ccd5..fd2d69c 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -37,9 +37,9 @@ The rest of the links [`ExpLink`](@ref), [`LogisticLink`](@ref), etc., are aliases for the corresponding wrapped functions in a `Link`. For example `ExpLink == Link{typeof(exp)}`. -When passing a [`Link`](@ref) to an `AbstractLikelihood`, this link +When passing a [`Link`](@ref) to a likelihood object, this link corresponds to the transformation `p=link(f)` while, as mentioned in the -[`Constrained parameters`](@ref) section, the statistics literature usually use +[Constrained parameters](@ref) section, the statistics literature usually uses the denomination [**inverse link or mean function**](https://en.wikipedia.org/wiki/Generalized_linear_model#Link_function) for it. ```@docs diff --git a/docs/src/index.md b/docs/src/index.md index 3b20c05..e1a2aad 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,6 +1,12 @@ ```@meta CurrentModule = GPLikelihoods ``` +```@setup test-repl +using Distributions +using GPLikelihoods +using StatsFuns +``` + # GPLikelihoods [`GPLikelihoods.jl`](https://github.com/JuliaGaussianProcesses/GPLikelihoods.jl) provides a practical interface to connect Gaussian and non-conjugate likelihoods @@ -14,11 +20,11 @@ taking a `Real` or an `AbstractVector` as an input and returns a Most likelihoods, like the [`GaussianLikelihood`](@ref), only require one latent Gaussian process. Passing a `Real` will therefore return a [`UnivariateDistribution`](https://juliastats.org/Distributions.jl/latest/univariate/), and passing an `AbstractVector{<:Real}` will return a [multivariate product of distributions](https://juliastats.org/Distributions.jl/latest/multivariate/#Product-distributions). -```@repl +```@repl test-repl f = 2.0; -GaussianLikelihood()(f) == Normal(2.0) -fs = [2.0, 3.0, 1.5] -GaussianLikelihood()(fs) == Product([Normal(2.0), Normal(3.0), Normal(1.5)]) +GaussianLikelihood()(f) == Normal(2.0, 1e-3) +fs = [2.0, 3.0, 1.5]; +GaussianLikelihood()(fs) isa AbstractMvNormal ``` Some likelihoods, like the [`CategoricalLikelihood`](@ref), requires multiple latent Gaussian processes, @@ -26,10 +32,10 @@ and an `AbstractVector{<:Real}` needs to be passed. To obtain a product of distributions an `AbstractVector{<:AbstractVector{<:Real}}` has to be passed (we recommend using [`ColVecs` and `RowVecs` from KernelFunctions.jl](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/api/#Vector-Valued-Inputs) if you need to transform an `AbstractMatrix`). -```@repl +```@repl test-repl fs = [2.0, 3.0, 4.5]; CategoricalLikelihood()(fs) isa Categorical -Fs = [rand(3) for _ in 1:4] +Fs = [rand(3) for _ in 1:4]; CategoricalLikelihood()(Fs) isa Product{<:Any,<:Categorical} ``` @@ -38,18 +44,18 @@ CategoricalLikelihood()(Fs) isa Product{<:Any,<:Categorical} The domain of some distributions parameters can be different from ``\mathbb{R}``, the real domain. To solve this problem, we also provide the [`Link`](@ref) type, which can be -passed to the [`Likelihood`](@ref) constructors. +passed to the likelihood constructors. Alternatively, `function`s can also directly be passed and will be wrapped in a `Link`). For more details about which likelihoods require a [`Link`](@ref) check out their docs. We typically named this passed link as the `invlink`. This comes from the statistic literature, where the "link" is defined as `f = link(y)`. -A classical example is the [`BernoulliLikelihood`](@ref) for classification, with the probability parameter ``p \in \[0, 1\]``. +A classical example is the [`BernoulliLikelihood`](@ref) for classification, with the probability parameter ``p \in [0, 1]``. The default it to use a [`logistic`](https://en.wikipedia.org/wiki/Logistic_function) transformation, but one could also use the inverse of the [`probit`](https://en.wikipedia.org/wiki/Probit) link: -```@repl +```@repl test-repl f = 2.0; BernoulliLikelihood()(f) == Bernoulli(logistic(f)) -BernoulliLikelihood(NormalCDFLink()) == Bernoulli(normalcdf(f)) +BernoulliLikelihood(NormalCDFLink())(f) == Bernoulli(normcdf(f)) ``` -Note that we passed the `inverse` of the `probit` function which is the `normalcdf` function. +Note that we passed the `inverse` of the `probit` function which is the `normcdf` function. diff --git a/src/links.jl b/src/links.jl index e7e66c9..515b8c1 100644 --- a/src/links.jl +++ b/src/links.jl @@ -8,6 +8,11 @@ A series of definitions are given in http://web.pdx.edu/~newsomj/mvclass/ho_link """ abstract type AbstractLink end +""" + ChainLink(links) + +Create a composed chain of different links. +""" struct ChainLink{Tls} <: AbstractLink links::Tls end From 0fc5a3632a52e539e4ee8718eb7060db00c3df0a Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Sat, 30 Apr 2022 01:41:56 +0000 Subject: [PATCH 40/44] CompatHelper: bump compat for "StatsFuns" to "1" --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 408e4ce..a198294 100644 --- a/Project.toml +++ b/Project.toml @@ -24,5 +24,5 @@ Functors = "0.1, 0.2" InverseFunctions = "0.1.2" IrrationalConstants = "0.1" SpecialFunctions = "1, 2" -StatsFuns = "0.9.13" +StatsFuns = "0.9.13, 1" julia = "1.6" From 7e5c4125f9800d7e87ea12223391ddc3ad9eca09 Mon Sep 17 00:00:00 2001 From: st-- Date: Mon, 2 May 2022 09:16:08 +0200 Subject: [PATCH 41/44] strict doc checks (#86) * turn on checkdocs * add missing docstrings * remove `elbo` reference * remove unknown keyword * fix footnote * light editing * clean up links explanation --- docs/make.jl | 8 +++++--- docs/src/api.md | 7 +++++++ docs/src/index.md | 29 ++++++++++++++++++----------- src/expectations.jl | 6 +++--- src/likelihoods/negativebinomial.jl | 4 ++-- 5 files changed, 35 insertions(+), 19 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index ce790be..dab5dc0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,13 +1,15 @@ using Documenter, GPLikelihoods makedocs(; - modules=[GPLikelihoods], + sitename="GPLikelihoods.jl", format=Documenter.HTML(), + modules=[GPLikelihoods], pages=["Home" => "index.md", "API" => "api.md"], repo="https://github.com/JuliaGaussianProcesses/GPLikelihoods.jl/blob/{commit}{path}#L{line}", - sitename="GPLikelihoods.jl", authors="JuliaGaussianProcesses organization", - assets=String[], + strict=true, + checkdocs=:exports, + #doctestfilters=JuliaGPsDocs.DOCTEST_FILTERS, ) deploydocs(; repo="github.com/JuliaGaussianProcesses/GPLikelihoods.jl", push_preview=true) diff --git a/docs/src/api.md b/docs/src/api.md index fd2d69c..79e4335 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -31,6 +31,7 @@ NBParamPower ```@docs Link ChainLink +BijectiveSimplexLink ``` The rest of the links [`ExpLink`](@ref), [`LogisticLink`](@ref), etc., @@ -54,3 +55,9 @@ ProbitLink NormalCDFLink SoftMaxLink ``` + +## Expectations + +```@docs +expected_loglikelihood +``` diff --git a/docs/src/index.md b/docs/src/index.md index e1a2aad..0c9c4c2 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -12,7 +12,7 @@ using StatsFuns [`GPLikelihoods.jl`](https://github.com/JuliaGaussianProcesses/GPLikelihoods.jl) provides a practical interface to connect Gaussian and non-conjugate likelihoods to Gaussian Processes. The API is very basic: Every `AbstractLikelihood` object is a [functor](https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects-1) -taking a `Real` or an `AbstractVector` as an input and returns a +taking a `Real` or an `AbstractVector` as an input and returning a `Distribution` from [`Distributions.jl`](https://github.com/JuliaStats/Distributions.jl). ### Single-latent vs multi-latent likelihoods @@ -27,7 +27,7 @@ fs = [2.0, 3.0, 1.5]; GaussianLikelihood()(fs) isa AbstractMvNormal ``` -Some likelihoods, like the [`CategoricalLikelihood`](@ref), requires multiple latent Gaussian processes, +Some likelihoods, like the [`CategoricalLikelihood`](@ref), require multiple latent Gaussian processes, and an `AbstractVector{<:Real}` needs to be passed. To obtain a product of distributions an `AbstractVector{<:AbstractVector{<:Real}}` has to be passed (we recommend using [`ColVecs` and `RowVecs` from KernelFunctions.jl](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/api/#Vector-Valued-Inputs) @@ -41,21 +41,28 @@ CategoricalLikelihood()(Fs) isa Product{<:Any,<:Categorical} ### Constrained parameters -The domain of some distributions parameters can be different from -``\mathbb{R}``, the real domain. -To solve this problem, we also provide the [`Link`](@ref) type, which can be -passed to the likelihood constructors. -Alternatively, `function`s can also directly be passed and will be wrapped in a `Link`). +The function values `f` of the latent Gaussian process live in the real domain +``\mathbb{R}``. For some likelihoods, the domain of the distribution parameter +`p` that is modulated by the latent Gaussian process is constrained to some +subset of ``\mathbb{R}``, e.g. only positive values or values in an interval. + +To connect these two domains, a transformation from `f` to `p` is required. +For this, we provide the [`Link`](@ref) type, which can be passed to the +likelihood constructors. (Alternatively, `function`s can also directly be +passed and will be wrapped in a `Link`.) + +We typically call this passed transformation the `invlink`. This comes from +the statistics literature, where the "link" is defined as `f = link(p)`, +whereas here we need `p = invlink(f)`. + For more details about which likelihoods require a [`Link`](@ref) check out their docs. -We typically named this passed link as the `invlink`. -This comes from the statistic literature, where the "link" is defined as `f = link(y)`. A classical example is the [`BernoulliLikelihood`](@ref) for classification, with the probability parameter ``p \in [0, 1]``. -The default it to use a [`logistic`](https://en.wikipedia.org/wiki/Logistic_function) transformation, but one could also use the inverse of the [`probit`](https://en.wikipedia.org/wiki/Probit) link: +The default is to use a [`logistic`](https://en.wikipedia.org/wiki/Logistic_function) transformation, but one could also use the inverse of the [`probit`](https://en.wikipedia.org/wiki/Probit) link: ```@repl test-repl f = 2.0; BernoulliLikelihood()(f) == Bernoulli(logistic(f)) BernoulliLikelihood(NormalCDFLink())(f) == Bernoulli(normcdf(f)) ``` -Note that we passed the `inverse` of the `probit` function which is the `normcdf` function. +Note that we passed the _inverse_ of the `probit` function which is the `normcdf` function. diff --git a/src/expectations.jl b/src/expectations.jl index 8855cfb..3de986b 100644 --- a/src/expectations.jl +++ b/src/expectations.jl @@ -42,11 +42,11 @@ callable that takes `f` as input and returns a Distribution over `y` that suppor ```math q(f) = ∫ p(f | u) q(u) du ``` -where `q(u)` is the variational distribution over inducing points (see -[`elbo`](@ref)). The marginal distributions of `q(f)` are given by `q_f`. +where `q(u)` is the variational distribution over inducing points. +The marginal distributions of `q(f)` are given by `q_f`. `quadrature` determines which method is used to calculate the expected log -likelihood - see [`elbo`](@ref) for more details. +likelihood. # Extended help diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl index 685528e..1c2ea40 100644 --- a/src/likelihoods/negativebinomial.jl +++ b/src/likelihoods/negativebinomial.jl @@ -9,7 +9,7 @@ default_invlink(::NBParamMean) = exp NegativeBinomialLikelihood(param::NBParam, invlink::Union{Function,Link}) There are many possible parametrizations for the Negative Binomial likelihood. -We follow the convention laid out in p.137 of [^1] and provide some common parametrizations. +We follow the convention laid out in p.137 of [^Hilbe'11] and provide some common parametrizations. The `NegativeBinomialLikelihood` has a special structure; the first type parameter `NBParam` defines in what parametrization the latent function is used, and contains the other (scalar) parameters. `NBParam` itself has two subtypes: @@ -51,7 +51,7 @@ julia> var(d) ≈ exp(2.0) * (1 + 3.0) true ``` -[^1] Hilbe, Joseph M. Negative binomial regression. Cambridge University Press, 2011. +[^Hilbe'11]: Hilbe, Joseph M. Negative binomial regression. Cambridge University Press, 2011. """ struct NegativeBinomialLikelihood{Tp<:NBParam,Tl<:AbstractLink} <: AbstractLikelihood params::Tp # Likelihood parametrization (and parameters) From bb12e890d1925828ab58b24dd9018d3fd732a32d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 10 May 2022 10:11:50 +0200 Subject: [PATCH 42/44] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a198294..ffc57b2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GPLikelihoods" uuid = "6031954c-0455-49d7-b3b9-3e1c99afaf40" authors = ["JuliaGaussianProcesses Team"] -version = "0.4.2" +version = "0.4.3" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From 38855fe68a530368caefd74e0e9e0e038b99c965 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 10 May 2022 10:12:23 +0200 Subject: [PATCH 43/44] Update Project.toml --- test/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Project.toml b/test/Project.toml index 9df0ed1..31bbe9e 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -9,6 +9,6 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] Distributions = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25" Functors = "0.1, 0.2" -StatsFuns = "0.9" +StatsFuns = "0.9, 1" Zygote = "0.6" julia = "1.3" From 6ad179dddf374fb4fa257be727bfe331ddf67943 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9o=20Galy-Fajou?= Date: Tue, 10 May 2022 10:12:46 +0200 Subject: [PATCH 44/44] Update Project.toml --- docs/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index bdd7e95..d6bd23b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,4 +8,4 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" Distributions = "0.25" Documenter = "0.25, 0.26, 0.27" GPLikelihoods = "0.4" -StatsFuns = "0.9" +StatsFuns = "0.9, 1"