Merge pull request #14 from JuliaGaussianProcesses/create_link

Adding a Link object

.JuliaFormatter.toml

@@ -0,0 +1 @@
+style = "blue"

.gitignore

@@ -3,6 +3,8 @@
 *.jl.mem
 .DS_Store
 /Manifest.toml
+test/Manifest.toml
 /dev/
 /docs/build/
 /docs/site/
+.vscode/

src/GPLikelihoods.jl

@@ -1,20 +1,33 @@
 module GPLikelihoods
 
-using Distributions
 using AbstractGPs
-using Random
+using Distributions
 using Functors
-using StatsFuns: logistic, softmax
-
-import Distributions
+using Random
+using StatsFuns
 
 export BernoulliLikelihood,
     CategoricalLikelihood,
-    GaussianLikelihood, 
-    HeteroscedasticGaussianLikelihood, 
+    GaussianLikelihood,
+    HeteroscedasticGaussianLikelihood,
     PoissonLikelihood,
     ExponentialLikelihood,
     GammaLikelihood
+export Link,
+    ChainLink,
+    ExpLink,
+    LogLink,
+    InvLink,
+    SqrtLink,
+    SquareLink,
+    LogitLink,
+    LogisticLink,
+    ProbitLink,
+    NormalCDFLink,
+    SoftMaxLink
+
+# Links
+include("links.jl")
 
 # Likelihoods
 include(joinpath("likelihoods", "bernoulli.jl"))

src/likelihoods/bernoulli.jl

@@ -1,16 +1,21 @@
 """
-    BernoulliLikelihood
+    BernoulliLikelihood(l::AbstractLink=LogisticLink())
 
 Bernoulli likelihood is to be used if we assume that the 
 uncertainity associated with the data follows a Bernoulli distribution.
+The link `l` needs to transform the input `f` to the domain [0, 1]
 
 ```math
-    p(y|f) = Bernoulli(y | f)
+    p(y|f) = Bernoulli(y | l(f))
 ```
-On calling, this would return a Bernoulli distribution with `f` probability of `true`.
+On calling, this would return a Bernoulli distribution with `l(f)` probability of `true`.
 """
-struct BernoulliLikelihood end
+struct BernoulliLikelihood{Tl<:AbstractLink}
+    invlink::Tl
+end
 
-(l::BernoulliLikelihood)(f::Real) = Bernoulli(logistic(f))
+BernoulliLikelihood() = BernoulliLikelihood(LogisticLink())
 
-(l::BernoulliLikelihood)(fs::AbstractVector{<:Real}) = Product(Bernoulli.(logistic.(fs)))
+(l::BernoulliLikelihood)(f::Real) = Bernoulli(l.invlink(f))
+
+(l::BernoulliLikelihood)(fs::AbstractVector{<:Real}) = Product(Bernoulli.(l.invlink.(fs)))

src/likelihoods/categorical.jl

@@ -1,16 +1,22 @@
 """
-    CategoricalLikelihood
+    CategoricalLikelihood(l::AbstractLink=SoftMaxLink())
 
 Categorical likelihood is to be used if we assume that the 
 uncertainity associated with the data follows a Categorical distribution.
 ```math
-    p(y|f_1, f_2, \\dots, f_{n-1}) = Categorical(y | softmax(f_1, f_2, \\dots, f_{n-1}, 0))
+    p(y|f_1, f_2, \\dots, f_{n-1}) = Categorical(y | l(f_1, f_2, \\dots, f_{n-1}, 0))
 ```
-On calling, this would return a Categorical distribution with `f_i` 
-probability of `i` category.
+Given an `AbstractVector` [f_1, f_2, ..., f_{n-1}], returns a `Categorical` distribution,
+with probabilities given by `l(f_1, f_2, ..., f_{n-1}, 0)`.
 """
-struct CategoricalLikelihood end
+struct CategoricalLikelihood{Tl<:AbstractLink}
+    invlink::Tl
+end
 
-(l::CategoricalLikelihood)(f::AbstractVector{<:Real}) = Categorical(softmax(vcat(f, 0)))
+CategoricalLikelihood() = CategoricalLikelihood(SoftMaxLink())
 
-(l::CategoricalLikelihood)(fs::AbstractVector) = Product(Categorical.(softmax.(vcat.(fs, 0))))
+(l::CategoricalLikelihood)(f::AbstractVector{<:Real}) = Categorical(l.invlink(vcat(f, 0)))
+
+function (l::CategoricalLikelihood)(fs::AbstractVector)
+    return Product(Categorical.(l.invlink.(vcat.(fs, 0))))
+end

src/likelihoods/exponential.jl

@@ -1,14 +1,20 @@
 """
-    ExponentialLikelihood()
+    ExponentialLikelihood(l::AbstractLink=ExpLink())
 
-Exponential likelihood with scale given by `exp(f)`.
+Exponential likelihood with scale given by `l(f)`.
 
 ```math
-    p(y|f) = Exponential(y | exp(f))
+    p(y|f) = Exponential(y | l(f))
 ```
 """
-struct ExponentialLikelihood end
+struct ExponentialLikelihood{Tl<:AbstractLink}
+    invlink::Tl
+end
 
-(l::ExponentialLikelihood)(f::Real) = Exponential(exp(f))
+ExponentialLikelihood() = ExponentialLikelihood(ExpLink())
 
-(l::ExponentialLikelihood)(fs::AbstractVector{<:Real}) = Product(Exponential.(exp.(fs)))
+(l::ExponentialLikelihood)(f::Real) = Exponential(l.invlink(f))
+
+function (l::ExponentialLikelihood)(fs::AbstractVector{<:Real})
+    return Product(Exponential.(l.invlink.(fs)))
+end

src/likelihoods/gamma.jl

@@ -1,21 +1,24 @@
 """
-    GammaLikelihood(α)
+    GammaLikelihood(α::Real=1.0, l::AbstractLink=ExpLink())
 
 Gamma likelihood with fixed shape `α`.
 
 ```math
-    p(y|f) = Gamma(y | α, θ=exp(f))
+    p(y|f) = Gamma(y | α, l(f))
 ```
-On calling, this would return a gamma distribution with shape `α` and scale `exp(f)`.
+On calling, this would return a gamma distribution with shape `α` and scale `l(f)`.
 """
-struct GammaLikelihood{T<:Real}
+struct GammaLikelihood{T<:Real,Tl<:AbstractLink}
     α::T    # shape parameter
+    invlink::Tl
 end
 
-GammaLikelihood() = GammaLikelihood(1.)
+GammaLikelihood() = GammaLikelihood(1.0)
+
+GammaLikelihood(α::Real) = GammaLikelihood(α, ExpLink())
 
 @functor GammaLikelihood
 
-(l::GammaLikelihood)(f::Real) = Gamma(l.α, exp(f))
+(l::GammaLikelihood)(f::Real) = Gamma(l.α, l.invlink(f))
 
-(l::GammaLikelihood)(fs::AbstractVector{<:Real}) = Product(Gamma.(l.α, exp.(fs)))
+(l::GammaLikelihood)(fs::AbstractVector{<:Real}) = Product(Gamma.(l.α, l.invlink.(fs)))

src/likelihoods/gaussian.jl

@@ -1,7 +1,7 @@
 """
     GaussianLikelihood(σ²)
 
-Gaussian likelihood with `σ²` variance. This is to be used if we assume that the 
+Gaussian likelihood with `σ²` variance. This is to be used if we assume that the
 uncertainity associated with the data follows a Gaussian distribution.
 
 ```math
@@ -10,31 +10,42 @@ uncertainity associated with the data follows a Gaussian distribution.
 On calling, this would return a normal distribution with mean `f` and variance σ².
 """
 struct GaussianLikelihood{T<:Real}
-    σ²::T
+    σ²::Vector{T}
 end
 
 GaussianLikelihood() = GaussianLikelihood(1e-6)
 
+GaussianLikelihood(σ²::Real) = GaussianLikelihood([σ²])
+
 @functor GaussianLikelihood
 
-(l::GaussianLikelihood)(f::Real) = Normal(f, sqrt(l.σ²))
+(l::GaussianLikelihood)(f::Real) = Normal(f, sqrt(first(l.σ²)))
 
-(l::GaussianLikelihood)(fs::AbstractVector{<:Real}) = MvNormal(fs, sqrt(l.σ²))
+(l::GaussianLikelihood)(fs::AbstractVector{<:Real}) = MvNormal(fs, sqrt(first(l.σ²)))
 
 """
-    HeteroscedasticGaussianLikelihood(σ²)
+    HeteroscedasticGaussianLikelihood(l::AbstractLink=ExpLink())
 
 Heteroscedastic Gaussian likelihood. 
 This is a Gaussian likelihood whose mean and the log of whose variance are functions of the
 latent process.
 
 ```math
-    p(y|[f, g]) = Normal(y | f, exp(g))
+    p(y|[f, g]) = Normal(y | f, l(g))
 ```
-On calling, this would return a normal distribution with mean `f` and variance `exp(g)`.
+On calling, this would return a normal distribution with mean `f` and std. dev. `l(g)`.
+Where `l` is link going from R to R^+
 """
-struct HeteroscedasticGaussianLikelihood end
+struct HeteroscedasticGaussianLikelihood{Tl<:AbstractLink}
+    invlink::Tl
+end
 
-(::HeteroscedasticGaussianLikelihood)(f::AbstractVector{<:Real}) = Normal(f[1], exp(f[2]))
+HeteroscedasticGaussianLikelihood() = HeteroscedasticGaussianLikelihood(ExpLink())
 
-(::HeteroscedasticGaussianLikelihood)(fs::AbstractVector) = MvNormal(first.(fs), exp.(last.(fs)))
+function (l::HeteroscedasticGaussianLikelihood)(f::AbstractVector{<:Real})
+    return Normal(f[1], l.invlink(f[2]))
+end
+
+function (l::HeteroscedasticGaussianLikelihood)(fs::AbstractVector)
+    return MvNormal(first.(fs), l.invlink.(last.(fs)))
+end

src/likelihoods/poisson.jl

@@ -1,11 +1,21 @@
 """
-    PoissonLikelihood()
+    PoissonLikelihood(l::AbstractLink=ExpLink())
 
-Poisson likelihood with rate as exponential of samples from GP `f`. This is to be used if
-we assume that the uncertainity associated with the data follows a Poisson distribution.
+Poisson likelihood with rate defined as `l(f)`.
+
+```math
+    p(y|f) = Poisson(y | θ=l(f))
+```
+
+This is to be used if  we assume that the uncertainity associated
+with the data follows a Poisson distribution.
 """
-struct PoissonLikelihood end
+struct PoissonLikelihood{L<:AbstractLink}
+    invlink::L
+end
+
+PoissonLikelihood() = PoissonLikelihood(ExpLink())
 
-(l::PoissonLikelihood)(f::Real) = Poisson(exp(f))
+(l::PoissonLikelihood)(f::Real) = Poisson(l.invlink(f))
 
-(l::PoissonLikelihood)(fs::AbstractVector{<:Real}) = Product(Poisson.(exp.(fs)))
+(l::PoissonLikelihood)(fs::AbstractVector{<:Real}) = Product(Poisson.(l.invlink.(fs)))

src/links.jl

@@ -0,0 +1,140 @@
+"""
+    AbstractLink
+
+Abstract type defining maps from R^n -> X.
+They can be applied by calling `link(x)`.
+
+A series of definitions are given in http://web.pdx.edu/~newsomj/mvclass/ho_link.pdf
+"""
+abstract type AbstractLink end
+
+struct ChainLink{Tls} <: AbstractLink
+    links::Tls
+end
+
+(l::ChainLink)(x) = foldl((x, l) -> l(x), l.ls; init=x)
+
+"""
+    Link(f)
+
+General construction for a link with a function `f`.
+"""
+struct Link{F} <: AbstractLink
+    f::F
+end
+
+(l::Link)(x) = l.f(x)
+
+"""
+    LogLink()
+
+`log` link, f:ℝ⁺->ℝ . Its inverse is the [`ExpLink`](@ref).
+"""
+struct LogLink <: AbstractLink end
+
+(::LogLink)(x) = log(x)
+
+Base.inv(::LogLink) = ExpLink()
+
+"""
+    ExpLink()
+
+`exp` link, f:ℝ->ℝ⁺. Its inverse is the [`LogLink`](@ref).
+"""
+struct ExpLink <: AbstractLink end
+
+(::ExpLink)(x) = exp(x)
+
+Base.inv(::ExpLink) = LogLink()
+
+"""
+    InvLink()
+
+`inv` link, f:ℝ/{0}->ℝ/{0}. It is its own inverse.
+"""
+struct InvLink <: AbstractLink end
+
+(::InvLink)(x) = inv(x)
+
+Base.inv(::InvLink) = InvLink()
+
+"""
+    SqrtLink()
+
+`sqrt` link, f:ℝ⁺∪{0}->ℝ⁺∪{0}. Its inverse is the [`SquareLink`](@ref).
+"""
+struct SqrtLink <: AbstractLink end
+
+(::SqrtLink)(x) = sqrt(x)
+
+Base.inv(::SqrtLink) = SquareLink()
+
+"""
+    SquareLink()
+
+`^2` link, f:ℝ->ℝ⁺∪{0}. Its inverse is the [`SqrtLink`](@ref).
+"""
+struct SquareLink <: AbstractLink end
+
+(::SquareLink)(x) = x^2
+
+Base.inv(::SquareLink) = SqrtLink()
+
+"""
+    LogitLink()
+
+`log(x/(1-x))` link, f:[0,1]->ℝ. Its inverse is the [`LogisticLink`](@ref).
+"""
+struct LogitLink <: AbstractLink end
+
+(::LogitLink)(x) = logit(x)
+
+Base.inv(::LogitLink) = LogisticLink()
+
+"""
+    LogisticLink()
+
+`exp(x)/(1+exp(-x))` link. f:ℝ->[0,1]. Its inverse is the [`Logit`](@ref).
+"""
+struct LogisticLink <: AbstractLink end
+
+(::LogisticLink)(x) = logistic(x)
+
+Base.inv(::LogisticLink) = LogitLink()
+
+"""
+    ProbitLink()
+
+`ϕ⁻¹(y)` link, where `ϕ⁻¹` is the `invcdf` of a `Normal` distribution, f:[0,1]->ℝ.
+Its inverse is the [`NormalCDFLink`](@ref).
+"""
+struct ProbitLink <: AbstractLink end
+
+(::ProbitLink)(x) = norminvcdf(x)
+
+Base.inv(::ProbitLink) = NormalCDFLink()
+
+"""
+    NormalCDFLink()
+
+`ϕ(y)` link, where `ϕ` is the `cdf` of a `Normal` distribution, f:ℝ->[0,1].
+Its inverse is the [`ProbitLink`](@ref).
+"""
+struct NormalCDFLink <: AbstractLink end
+
+(::NormalCDFLink)(x) = normcdf(x)
+
+Base.inv(::NormalCDFLink) = ProbitLink()
+
+(::ChainLink{<:Tuple{LogLink,NormalCDFLink}})(x) = normlogcdf(x) # Specialisation for log + normal cdf
+
+"""
+    SoftMaxLink()
+
+`softmax` link, i.e `f(x)ᵢ = exp(xᵢ)/∑ⱼexp(xⱼ)`.
+f:ℝⁿ->Sⁿ⁻¹, where Sⁿ⁻¹ is an [(n-1)-simplex](https://en.wikipedia.org/wiki/Simplex)
+It has no defined inverse
+"""
+struct SoftMaxLink <: AbstractLink end
+
+(::SoftMaxLink)(x::AbstractVector{<:Real}) = softmax(x)