Add BernoulliLogit

docs/src/univariate.md

@@ -473,6 +473,7 @@ plotdensity((0.001, 3), Weibull, (0.5, 1)) # hide
 
 ```@docs
 Bernoulli
+BernoulliLogit
 BetaBinomial
 Binomial
 Categorical

src/Distributions.jl

@@ -67,6 +67,7 @@ export
     # distribution types
     Arcsine,
     Bernoulli,
+    BernoulliLogit,
     Beta,
     BetaBinomial,
     BetaPrime,

src/univariate/discrete/bernoullilogit.jl

@@ -0,0 +1,108 @@
+"""
+    BernoulliLogit(logitp=0.0)
+
+A *Bernoulli distribution* that is parameterized by the logit `logitp = logit(p) = log(p/(1-p))` of its success rate `p`.
+
+```math
+P(X = k) = \\begin{cases}
+\\operatorname{logistic}(-logitp) = \\frac{1}{1 + \\exp{(logitp)}} & \\quad \\text{for } k = 0, \\\\
+\\operatorname{logistic}(logitp) = \\frac{1}{1 + \\exp{(-logitp)}} & \\quad \\text{for } k = 1.
+\\end{cases}
+```
+
+External links:
+
+* [Bernoulli distribution on Wikipedia](http://en.wikipedia.org/wiki/Bernoulli_distribution)
+
+See also [`Bernoulli`](@ref)
+"""
+struct BernoulliLogit{T<:Real} <: DiscreteUnivariateDistribution
+    logitp::T
+end
+
+BernoulliLogit() = BernoulliLogit(0.0)
+
+@distr_support BernoulliLogit false true
+
+Base.eltype(::Type{<:BernoulliLogit}) = Bool
+
+#### Conversions
+Base.convert(::Type{BernoulliLogit{T}}, d::BernoulliLogit) where {T<:Real} = BernoulliLogit{T}(T(d.logitp))
+Base.convert(::Type{BernoulliLogit{T}}, d::BernoulliLogit{T}) where {T<:Real} = d
+
+#### Parameters
+
+succprob(d::BernoulliLogit) = logistic(d.logitp)
+failprob(d::BernoulliLogit) = logistic(-d.logitp)
+logsuccprob(d::BernoulliLogit) = -log1pexp(-d.logitp)
+logfailprob(d::BernoulliLogit) = -log1pexp(d.logitp)
+
+params(d::BernoulliLogit) = (d.logitp,)
+partype(::BernoulliLogit{T}) where {T} = T
+
+#### Properties
+
+mean(d::BernoulliLogit) = succprob(d)
+var(d::BernoulliLogit) =  succprob(d) * failprob(d)
+function skewness(d::BernoulliLogit)
+    p0 = failprob(d)
+    p1 = succprob(d)
+    return (p0 - p1) / sqrt(p0 * p1)
+end
+kurtosis(d::BernoulliLogit) = 1 / var(d) - 6
+
+mode(d::BernoulliLogit) = d.logitp > 0 ? 1 : 0
+
+function modes(d::BernoulliLogit)
+    logitp = d.logitp
+    z = zero(logitp)
+    logitp < z ? [false] : (logitp > z ? [true] : [false, true])
+end
+
+median(d::BernoulliLogit) = d.logitp > 0
+
+function entropy(d::BernoulliLogit)
+    logitp = d.logitp
+    (logitp == -Inf || logitp == Inf) ? float(zero(logitp)) : (logitp > 0 ? -(succprob(d) * logitp + logfailprob(d)) : -(logsuccprob(d) - failprob(d) * logitp))
+end
+
+#### Evaluation
+
+pdf(d::BernoulliLogit, x::Bool) = x ? succprob(d) : failprob(d)
+pdf(d::BernoulliLogit, x::Real) = x == 0 ? failprob(d) : (x == 1 ? succprob(d) : zero(float(d.logitp)))
+
+logpdf(d::BernoulliLogit, x::Bool) = x ? logsuccprob(d) : logfailprob(d)
+logpdf(d::BernoulliLogit, x::Real) = x == 0 ? logfailprob(d) : (x == 1 ? logsuccprob(d) : oftype(float(d.logitp), -Inf))
+
+cdf(d::BernoulliLogit, x::Bool) = x ? one(float(d.logitp)) : failprob(d)
+cdf(d::BernoulliLogit, x::Int) = x < 0 ? zero(float(d.logitp)) : (x < 1 ? failprob(d) : one(float(d.logitp)))
+
+logcdf(d::BernoulliLogit, x::Bool) = x ? zero(float(d.logitp)) : logfailprob(d)
+logcdf(d::BernoulliLogit, x::Int) = x < 0 ? oftype(float(d.logitp), -Inf) : (x < 1 ? logfailprob(d) : zero(float(d.logitp)))
+
+ccdf(d::BernoulliLogit, x::Bool) = x ? zero(float(d.logitp)) : succprob(d)
+ccdf(d::BernoulliLogit, x::Int) = x < 0 ? one(float(d.logitp)) : (x < 1 ? succprob(d) : zero(float(d.logitp)))
+
+logccdf(d::BernoulliLogit, x::Bool) = x ? oftype(float(d.logitp), -Inf) : logsuccprob(d)
+logccdf(d::BernoulliLogit, x::Int) = x < 0 ? zero(float(d.logitp)) : (x < 1 ? logsuccprob(d) : oftype(float(d.logitp), -Inf))
+
+function quantile(d::BernoulliLogit, p::Real)
+    T = float(partype(d))
+    0 <= p <= 1 ? (p <= failprob(d) ? zero(T) : one(T)) : T(NaN)
+end
+function cquantile(d::BernoulliLogit, p::Real)
+    T = float(partype(d))
+    0 <= p <= 1 ? (p >= succprob(d) ? zero(T) : one(T)) : T(NaN)
+end
+
+mgf(d::BernoulliLogit, t::Real) = failprob(d) + exp(t + logsuccprob(d))
+function cgf(d::BernoulliLogit, t)
+    # log(1-p+p*exp(t)) = logaddexp(log(1-p), t + log(p))
+    logaddexp(logfailprob(d), t + logsuccprob(d))
+end
+cf(d::BernoulliLogit, t::Real) = failprob(d) + succprob(d) * cis(t)
+
+
+#### Sampling
+
+rand(rng::AbstractRNG, d::BernoulliLogit) = logit(rand(rng)) <= d.logitp

src/univariates.jl

@@ -650,6 +650,7 @@ end
 
 const discrete_distributions = [
     "bernoulli",
+    "bernoullilogit",
     "betabinomial",
     "binomial",
     "dirac",

test/univariate/discrete/bernoullilogit.jl

@@ -0,0 +1,81 @@
+using Distributions
+using StatsFuns
+using Test, Random
+
+@testset "basic properties" begin
+    @test BernoulliLogit() === BernoulliLogit(0.0)
+
+    for logitp in (-0.3, 0.2, 0.1f0)
+        d = BernoulliLogit(logitp)
+        @test d isa BernoulliLogit{typeof(logitp)}
+        @test convert(typeof(d), d) === d
+        @test convert(BernoulliLogit{Float16}, d) === BernoulliLogit(Float16(logitp))
+        @test eltype(typeof(d)) === Bool
+        @test params(d) == (logitp,)
+        @test partype(d) === typeof(logitp)
+    end
+end
+
+@testset "succprob/failprob" begin
+    for p in (0.0, 0.1, 0.31f0, 0.5, 0.7f0, 0.95, 1.0)
+        d = BernoulliLogit(logit(p))
+        @test @inferred(succprob(d)) ≈ p
+        @test @inferred(failprob(d)) ≈ 1 - p
+        @test @inferred(Distributions.logsuccprob(d)) ≈ log(p)
+        @test @inferred(Distributions.logfailprob(d)) ≈ log1p(-p)
+    end
+end
+
+@testset "rand" begin
+    @test rand(BernoulliLogit()) isa Bool
+    @test rand(BernoulliLogit(), 10) isa Vector{Bool}
+
+    N = 10_000
+    for p in (0.0, 0.1, 0.31f0, 0.5, 0.7f0, 0.95, 1.0)
+        d = BernoulliLogit(logit(p))
+        @test @inferred(rand(d)) isa Bool
+        @test @inferred(rand(d, 10)) isa Vector{Bool}
+        @test mean(rand(d, N)) ≈ p atol=0.01
+    end
+end
+
+@testset "cgf" begin
+    test_cgf(BernoulliLogit(), (1f0, -1f0, 1e6, -1e6))
+    test_cgf(BernoulliLogit(0.1), (1f0, -1f0, 1e6, -1e6))
+end
+
+@testset "comparison with `Bernoulli`" begin
+    for p in (0.0, 0.1, 0.31f0, 0.5, 0.7f0, 0.95, 1.0)
+        d = BernoulliLogit(logit(p))
+        d0 = Bernoulli(p)
+
+        @test @inferred(mean(d)) ≈ mean(d0)
+        @test @inferred(var(d)) ≈ var(d0)
+        @test @inferred(skewness(d)) ≈ skewness(d0)
+        @test @inferred(kurtosis(d)) ≈ kurtosis(d0)
+        @test @inferred(mode(d)) ≈ mode(d0)
+        @test @inferred(modes(d)) ≈ modes(d0)
+        @test @inferred(median(d)) ≈ median(d0)
+        @test @inferred(entropy(d)) ≈ entropy(d0)
+
+        for x in (true, false, 0, 1, -3, 5)
+            @test @inferred(pdf(d, x)) ≈ pdf(d0, x)
+            @test @inferred(logpdf(d, x)) ≈ logpdf(d0, x)
+            @test @inferred(cdf(d, x)) ≈ cdf(d0, x)
+            @test @inferred(logcdf(d, x)) ≈ logcdf(d0, x)
+            @test @inferred(ccdf(d, x)) ≈ ccdf(d0, x)
+            @test @inferred(logccdf(d, x)) ≈ logccdf(d0, x)
+        end
+
+        for q in (-0.2f0, 0.25, 0.6f0, 1.5)
+            @test @inferred(quantile(d, q)) ≈ quantile(d0, q) nans=true
+            @test @inferred(cquantile(d, q)) ≈ cquantile(d0, q) nans=true
+        end
+
+        for t in (-5.2, 1.2f0)
+            @test @inferred(mgf(d, t)) ≈ mgf(d0, t) rtol=1e-6
+            @test @inferred(cgf(d, t)) ≈ cgf(d0, t) rtol=1e-6
+            @test @inferred(cf(d, t)) ≈ cf(d0, t) rtol=1e-6
+        end
+    end
+end