Merge branch 'master' into tgf/softmax_docs

Project.toml

@@ -1,7 +1,7 @@
 name = "GPLikelihoods"
 uuid = "6031954c-0455-49d7-b3b9-3e1c99afaf40"
 authors = ["JuliaGaussianProcesses Team"]
-version = "0.4.1"
+version = "0.4.3"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -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"

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, 1"

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)

docs/src/api.md

@@ -15,6 +15,17 @@ HeteroscedasticGaussianLikelihood
 PoissonLikelihood
 ```
 
+### Negative Binomial
+
+```@docs
+NegativeBinomialLikelihood
+NBParamSuccess
+NBParamFailure
+NBParamI
+NBParamII
+NBParamPower
+```
+
 ## Links
 
 ```@docs
@@ -27,9 +38,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
@@ -43,4 +54,10 @@ LogisticLink
 ProbitLink
 NormalCDFLink
 SoftMaxLink
-```
+```
+
+## Expectations
+
+```@docs
+expected_loglikelihood
+```

docs/src/index.md

@@ -1,55 +1,68 @@
 ```@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
 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
 
 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,
+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)
 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}
 ```
 
 ### 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`](@ref) 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:
+A classical example is the [`BernoulliLikelihood`](@ref) for classification, with the probability parameter ``p \in [0, 1]``.
+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
+```@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.

src/GPLikelihoods.jl

@@ -15,6 +15,7 @@ export BernoulliLikelihood,
     ExponentialLikelihood,
     GammaLikelihood,
     NegativeBinomialLikelihood
+export NBParamI, NBParamII, NBParamPower, NBParamSuccess, NBParamFailure
 export Link,
     ChainLink,
     BijectiveSimplexLink,

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
 

src/likelihoods/negativebinomial.jl

@@ -1,36 +1,169 @@
+abstract type NBParam end
+
+abstract type NBParamProb <: NBParam end
+default_invlink(::NBParamProb) = logistic
+abstract type NBParamMean <: NBParam end
+default_invlink(::NBParamMean) = exp
+
 """
-    NegativeBinomialLikelihood(l=logistic; successes::Real=1)
+    NegativeBinomialLikelihood(param::NBParam, invlink::Union{Function,Link})
 
-Negative binomial likelihood with number of successes `successes`.
+There are many possible parametrizations for the Negative Binomial likelihood.
+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:
+- `NBParamProb` for parametrizations where `f -> p`, the probability of success of a Bernoulli event
+- `NBParamMean` for parametrizations where `f -> μ`, the expected number of events
 
-```math
-    p(k|successes, f) = \\frac{\\Gamma(k+successes)}{k! \\Gamma(successes)} l(f)^successes (1 - l(f))^k
+## `NBParam` predefined types
+
+### `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 `μ = 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 `μ(1 + α * μ^ρ)`
+
+
+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
+
+```julia-repl
+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
 ```
-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
+
+[^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)
     invlink::Tl
+    function NegativeBinomialLikelihood(
+        params::Tparam, invlink=default_invlink(params)
+    ) where {Tparam<:NBParam}
+        # 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
 end
 
-function NegativeBinomialLikelihood(l=logistic; successes::Real=1)
-    return NegativeBinomialLikelihood(successes, link(l))
+@deprecate NegativeBinomialLikelihood(l=logistic; successes=1) NegativeBinomialLikelihood(
+    NBParamI(successes), l
+)
+
+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)`. ",
+        "For a full explanation, see `NegativeBinomialLikelihood` docs",
+    )
 end
 
 @functor NegativeBinomialLikelihood
 
-(l::NegativeBinomialLikelihood)(f::Real) = NegativeBinomial(l.successes, l.invlink(f))
-
 (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs))
+
+@doc raw"""
+    NBParamSuccess(successes)
+
+Negative Binomial parametrization with `successes` the number of successes and
+`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(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
+    successes::T
+end
+
+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
+`invlink(f)` the probability of success.
+This corresponds to the definition used by [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution).
+
+```math
+  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
+    failures::T
+end
+
+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_excessvar_to_r_p(μ::Real, ev::Real) = μ / ev, 1 / (1 + ev)
+
+"""
+    NBParamI(α)
+
+Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α)` where `α > 0`.
+"""
+struct NBParamI{T} <: NBParamMean
+    α::T
+end
+
+function (l::NegativeBinomialLikelihood{<:NBParamI})(f::Real)
+    μ = l.invlink(f)
+    ev = l.params.α
+    return NegativeBinomial(_nb_mean_excessvar_to_r_p(μ, ev)...)
+end
+
+"""
+    NBParamII(α)
+
+Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α * μ)` where `α > 0`.
+"""
+struct NBParamII{T} <: NBParamMean
+    α::T
+end
+
+function (l::NegativeBinomialLikelihood{<:NBParamII})(f::Real)
+    μ = l.invlink(f)
+    ev = l.params.α * μ
+    return NegativeBinomial(_nb_mean_excessvar_to_r_p(μ, ev)...)
+end
+
+"""
+    NBParamPower(α, ρ)
+
+Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α * μ^ρ)` where `α > 0` and `ρ > 0`.
+"""
+struct NBParamPower{Tα,Tρ} <: NBParamMean
+    α::Tα
+    ρ::Tρ
+end
+
+function (l::NegativeBinomialLikelihood{<:NBParamPower})(f::Real)
+    μ = l.invlink(f)
+    ev = l.params.α * μ^l.params.ρ
+    return NegativeBinomial(_nb_mean_excessvar_to_r_p(μ, ev)...)
+end