Add Weibull

CHANGELOG.md

@@ -15,7 +15,8 @@ and this project adheres to [Semantic Versioning].
 - The `RealLine`, `NonNegativeReals`, `PositiveReals`, `NonNegativeIntegers`,
   `RealInterval`, `IntegerRange`, `IntegerSimplex` and `RealVectors` supports.
   `UnitInterval` will arrive with the first measure that needs one.
-- `Normal(μ, σ)`, `LogNormal(μ, σ)`, `Exponential(θ)` and `Uniform(a, b)`, with
+- `Normal(μ, σ)`, `LogNormal(μ, σ)`, `Exponential(θ)`, `Weibull(α, θ)` and
+  `Uniform(a, b)`, with
   heterogeneous parameter types and no promotion or validation at construction.
 - `MvNormal(μ, L)`, the first multivariate measure, with the
   `ContinuousMultivariateMeasure` alias it dispatches on. It is parameterized by the

README.md

@@ -11,8 +11,8 @@ density and sampling operations, and compatible with automatic differentiation,
 broadcasting on GPU arrays, and Reactant tracing.
 
 The package is experimental. At present it implements `Normal`, `LogNormal`,
-`Exponential`, `Uniform`, `Laplace`, `Cauchy`, `Categorical`, `Bernoulli`, `Binomial`,
-`Poisson`, `Geometric`, `MvNormal`, and `Multinomial`.
+`Exponential`, `Weibull`, `Uniform`, `Laplace`, `Cauchy`, `Categorical`, `Bernoulli`,
+`Binomial`, `Poisson`, `Geometric`, `MvNormal`, and `Multinomial`.
 
 ## Installation
 
@@ -61,7 +61,7 @@ than throwing.
 
 ## Available API
 
-`Normal(μ, σ)`, `LogNormal(μ, σ)`, `Exponential(θ)`, `Uniform(a, b)`,
+`Normal(μ, σ)`, `LogNormal(μ, σ)`, `Exponential(θ)`, `Weibull(α, θ)`, `Uniform(a, b)`,
 `Laplace(μ, b)`, `Cauchy(μ, σ)`, `Categorical(p)`, `Bernoulli(p)`, `Binomial(n, p)`,
 `Poisson(λ)`, and `Geometric(p)` each support:
 
@@ -205,9 +205,9 @@ See the [contribution guide](docs/src/90-contributing.md) for contribution guide
 ## Current scope
 
 ProbabilityMeasures.jl currently contains `Normal`, `LogNormal`, `Exponential`,
-`Uniform`, `Cauchy`, `Laplace`, `Categorical`, `Bernoulli`, `Binomial`, `Poisson`,
-`Geometric`, `MvNormal`, and `Multinomial`. Transformed or composite measures and
-Distributions.jl interoperability are not implemented yet.
+`Weibull`, `Uniform`, `Cauchy`, `Laplace`, `Categorical`, `Bernoulli`, `Binomial`,
+`Poisson`, `Geometric`, `MvNormal`, and `Multinomial`. Transformed or composite measures
+and Distributions.jl interoperability are not implemented yet.
 
 ## Citation
 

libs/ProbabilityMeasuresTest/src/implementations.jl

@@ -23,6 +23,16 @@ _exactparams(::LogNormal) = LogNormal(0, 1)
 _invalids(::Exponential) = (Exponential(-1.0), Exponential(0.0), Exponential(Inf))
 _exactparams(::Exponential) = Exponential(1)
 
+@implements MeasureInterface{UNIVARIATE_OPTIONALS} Weibull [
+    Weibull(1.0, 1.0), Weibull(0.75, 2.5), Weibull(2.0f0, 0.5f0)
+]
+
+function _invalids(::Weibull)
+    return (Weibull(-1.0, 1.0), Weibull(0.0, 1.0), Weibull(1.0, 0.0), Weibull(Inf, 1.0))
+end
+
+_exactparams(::Weibull) = Weibull(2, 3)
+
 @implements MeasureInterface{UNIVARIATE_OPTIONALS} Uniform [
     Uniform(0.0, 1.0), Uniform(-1.0, 2.0), Uniform(0.0f0, 2.0f0)
 ]

src/ProbabilityMeasures.jl

@@ -27,6 +27,7 @@ include("core/interface.jl")
 include("univariate/continuous/normal.jl")
 include("univariate/continuous/lognormal.jl")
 include("univariate/continuous/exponential.jl")
+include("univariate/continuous/weibull.jl")
 include("univariate/continuous/uniform.jl")
 include("univariate/continuous/laplace.jl")
 include("univariate/continuous/cauchy.jl")
@@ -74,6 +75,7 @@ export mean, var, std, median, quantile, cov
 export Normal
 export LogNormal
 export Exponential
+export Weibull
 export Uniform
 export Laplace
 export Cauchy

src/univariate/continuous/weibull.jl

@@ -0,0 +1,110 @@
+"""
+    Weibull(α, θ)
+    Weibull()
+
+The Weibull measure on ``[0, \\infty)`` with shape `α` and scale `θ`. Its density is
+
+```math
+p(x) = \\frac{\\alpha}{\\theta}
+       \\left(\\frac{x}{\\theta}\\right)^{\\alpha - 1}
+       \\exp\\!\\left[-\\left(\\frac{x}{\\theta}\\right)^\\alpha\\right]
+```
+
+`Weibull()` creates the unit-shape, unit-scale measure using `Float64` values.
+
+The constructor does not check its arguments. Invalid parameters give a non-finite
+density. Use [`checkparams`](@ref) to check them when needed.
+"""
+struct Weibull{A<:Number,T<:Number} <: ContinuousUnivariateMeasure
+    α::A
+    θ::T
+end
+
+Weibull() = Weibull(1.0, 1.0)
+
+Base.eltype(::Type{Weibull{A,T}}) where {A,T} = float(promote_type(A, T))
+
+function checkparams(d::Weibull)
+    return isfinite(d.α) & isfinite(d.θ) & (d.α > zero(d.α)) & (d.θ > zero(d.θ))
+end
+
+support(::Weibull) = NonNegativeReals()
+
+@inline sval(d::Weibull, x::Number) = x / d.θ
+
+@inline function DensityInterface.logdensityof(d::Weibull, x::Number)
+    scaled, shape = promote(sval(d, x), d.α)
+    s = abs(float(scaled))
+    α, θ = oftype(s, shape), oftype(s, d.θ)
+    value = logt(α) - logt(θ) + (α - one(α)) * logt(s) - s^α
+    atzero = select(
+        α < one(α),
+        () -> oftype(s, Inf),
+        () -> select(α == one(α), () -> -logt(θ), () -> oftype(s, -Inf)),
+    )
+    return select(
+        x > zero(x),
+        () -> value,
+        () -> select(x == zero(x), () -> atzero, () -> oftype(s, -Inf)),
+    )
+end
+
+@inline Base.rand(rng::AbstractRNG, d::Weibull) = quantile(d, rand(rng, noisetype(d)))
+
+function Statistics.mean(d::Weibull)
+    α, θ = promote(d.α, d.θ)
+    α = float(α)
+    return oftype(α, θ) * exp(loggamma(one(α) + inv(α)))
+end
+
+function Statistics.var(d::Weibull)
+    α, θ = promote(d.α, d.θ)
+    α, θ = float(α), oftype(float(α), θ)
+    g₁ = loggamma(one(α) + inv(α))
+    g₂ = loggamma(one(α) + 2 / α)
+    return θ^2 * exp(2g₁) * expm1(g₂ - 2g₁)
+end
+
+function entropy(d::Weibull)
+    α, θ = promote(d.α, d.θ)
+    α, θ = float(α), oftype(float(α), θ)
+    γ = oftype(α, Base.MathConstants.eulergamma)
+    return γ * (one(α) - inv(α)) + logt(θ / α) + one(α)
+end
+
+function cdf(d::Weibull, x::Number)
+    s, α = promote(sval(d, x), d.α)
+    t = abs(float(s))^oftype(float(s), α)
+    return select(x >= zero(x), () -> -expm1(-t), () -> zero(t))
+end
+
+function ccdf(d::Weibull, x::Number)
+    s, α = promote(sval(d, x), d.α)
+    t = abs(float(s))^oftype(float(s), α)
+    return select(x >= zero(x), () -> exp(-t), () -> one(t))
+end
+
+function logcdf(d::Weibull, x::Number)
+    s, α = promote(sval(d, x), d.α)
+    t = abs(float(s))^oftype(float(s), α)
+    return select(x >= zero(x), () -> log1mexpt(-t), () -> oftype(t, -Inf))
+end
+
+function logccdf(d::Weibull, x::Number)
+    s, α = promote(sval(d, x), d.α)
+    t = abs(float(s))^oftype(float(s), α)
+    return select(x >= zero(x), () -> -t, () -> zero(t))
+end
+
+function Statistics.quantile(d::Weibull, p::Number)
+    probability, shape, scale = promote(p, d.α, d.θ)
+    q = float(probability)
+    α, θ = oftype(q, shape), oftype(q, scale)
+    value = θ * abs(-log1pt(-q))^inv(α)
+    valid = (q >= zero(q)) & (q <= one(q))
+    return select(valid, () -> value, () -> oftype(q, NaN))
+end
+
+function Base.show(io::IO, d::Weibull)
+    return print(io, "Weibull(α=", d.α, ", θ=", d.θ, ")")
+end

test/reactant/runtests.jl

@@ -18,6 +18,7 @@ using Test
         Uniform(-1.0, 2.0),
         Cauchy(-1.0, 2.0),
         Geometric(0.3),
+        Weibull(1.5, 2.0),
     )
         test_reactant(d, default_testpoints(d))
     end
@@ -47,5 +48,10 @@ using Test
         grad = p -> Enzyme.gradient(Enzyme.Reverse, loss, p)
         got = @jit grad(Reactant.ConcreteRNumber(0.5))
         @test Float64(got[1]) ≈ -4.0 rtol = 1e-10
+
+        loss = (a, s) -> logdensityof(Weibull(a, s), 1.0)
+        grad = (a, s) -> Enzyme.gradient(Enzyme.Reverse, loss, a, s)
+        got = @jit grad(Reactant.ConcreteRNumber(1.0), Reactant.ConcreteRNumber(1.0))
+        @test [Float64(got[1]), Float64(got[2])] ≈ [1.0, 0.0] rtol = 1e-10
     end
 end

test/test-weibull.jl

@@ -0,0 +1,69 @@
+using ProbabilityMeasures
+using ProbabilityMeasuresTest: test_measure
+using Distributions: Distributions
+using ForwardDiff: ForwardDiff
+using Random: Xoshiro
+using Test
+
+@testset "conformance" begin
+    reference_logpdf(m, x) =
+        Distributions.logpdf(Distributions.Weibull(Float64(m.α), Float64(m.θ)), x)
+    for d in (Weibull(1.0, 1.0), Weibull(0.75, 2.5), Weibull(2.0f0, 0.5f0), Weibull(2, 3))
+        test_measure(d; name=string(d), reference_logpdf=reference_logpdf)
+    end
+end
+
+@testset "parameters and support" begin
+    dual = ForwardDiff.Dual(2.0, 1.0)
+    @test typeof(Weibull(dual, 1.0)) === Weibull{typeof(dual),Float64}
+    @test typeof(Weibull(2.0f0, 1)) === Weibull{Float32,Int}
+    @test support(Weibull()) === NonNegativeReals()
+
+    @test checkparams(Weibull(2.0, 3.0))
+end
+
+@testset "density at zero" begin
+    @test logdensityof(Weibull(0.5, 2.0), 0.0) == Inf
+    @test logdensityof(Weibull(1.0, 2.0), 0.0) == -log(2.0)
+    @test logdensityof(Weibull(2.0, 2.0), 0.0) == -Inf
+    @test logdensityof(Weibull(0.5, 2.0), -1.0) == -Inf
+end
+
+@testset "distribution-function tails" begin
+    exponential = Weibull(1.0, 1.0)
+    @test cdf(exponential, 1e-20) ≈ 1e-20
+    @test ccdf(exponential, 1000.0) == 0.0
+    @test logccdf(exponential, 1000.0) == -1000.0
+
+    for p in (1e-20, prevfloat(1.0))
+        @test cdf(exponential, quantile(exponential, p)) ≈ p
+    end
+end
+
+@testset "reference numerics against Distributions.jl" begin
+    for (α, θ) in ((0.5, 0.4), (1.0, 1.0), (2.5, 3.0)), x in (0.2, 1.7, 8.0)
+        d, r = Weibull(α, θ), Distributions.Weibull(α, θ)
+        @test logdensityof(d, x) ≈ Distributions.logpdf(r, x)
+        @test cdf(d, x) ≈ Distributions.cdf(r, x)
+        @test ccdf(d, x) ≈ Distributions.ccdf(r, x)
+        @test logcdf(d, x) ≈ Distributions.logcdf(r, x)
+        @test logccdf(d, x) ≈ Distributions.logccdf(r, x)
+    end
+    for (α, θ) in ((0.5, 0.4), (1.0, 1.0), (2.5, 3.0)), p in (0.01, 0.5, 0.999)
+        @test quantile(Weibull(α, θ), p) ≈
+            Distributions.quantile(Distributions.Weibull(α, θ), p)
+    end
+    for (α, θ) in ((0.5, 0.4), (1.0, 1.0), (2.5, 3.0))
+        d, r = Weibull(α, θ), Distributions.Weibull(α, θ)
+        @test mean(d) ≈ Distributions.mean(r)
+        @test var(d) ≈ Distributions.var(r)
+        @test std(d) ≈ Distributions.std(r)
+        @test entropy(d) ≈ Distributions.entropy(r)
+    end
+end
+
+@testset "inverse-CDF sampling" begin
+    d = Weibull(1.5, 2.0)
+    seed = 0x4153554b41
+    @test rand(Xoshiro(seed), d) == quantile(d, rand(Xoshiro(seed), Float64))
+end