From c4c0f362182a521516fad95baecabb7c8ee91149 Mon Sep 17 00:00:00 2001 From: Ben Cox Date: Sat, 22 Aug 2026 22:52:14 +0200 Subject: [PATCH] Add Weibull --- CHANGELOG.md | 3 +- README.md | 12 +- .../src/implementations.jl | 10 ++ src/ProbabilityMeasures.jl | 2 + src/univariate/continuous/weibull.jl | 110 ++++++++++++++++++ test/reactant/runtests.jl | 6 + test/test-weibull.jl | 69 +++++++++++ 7 files changed, 205 insertions(+), 7 deletions(-) create mode 100644 src/univariate/continuous/weibull.jl create mode 100644 test/test-weibull.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 8d0d36e..003be80 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/README.md b/README.md index ec216e7..7f34a1a 100644 --- a/README.md +++ b/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 diff --git a/libs/ProbabilityMeasuresTest/src/implementations.jl b/libs/ProbabilityMeasuresTest/src/implementations.jl index a3ce857..5dc0dc6 100644 --- a/libs/ProbabilityMeasuresTest/src/implementations.jl +++ b/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) ] diff --git a/src/ProbabilityMeasures.jl b/src/ProbabilityMeasures.jl index df328a8..1956311 100644 --- a/src/ProbabilityMeasures.jl +++ b/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 diff --git a/src/univariate/continuous/weibull.jl b/src/univariate/continuous/weibull.jl new file mode 100644 index 0000000..cff8af5 --- /dev/null +++ b/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 diff --git a/test/reactant/runtests.jl b/test/reactant/runtests.jl index d36f8ff..6f8ac74 100644 --- a/test/reactant/runtests.jl +++ b/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 diff --git a/test/test-weibull.jl b/test/test-weibull.jl new file mode 100644 index 0000000..70a0338 --- /dev/null +++ b/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