Irwin-Hall and Bates distributions

docs/src/bestiary/continuous.md

@@ -111,6 +111,12 @@ Gompertz
 plotdensity((0.0, 1.0), Gompertz, (1,1)) # hide
 ```
 ```@docs
+IrwinHall
+```
+```@example plotdensity
+plotdensity((0.0, 3.0), IrwinHall, (3)) # hide
+```
+```@docs
 Lomax
 ```
 ```@example plotdensity

src/AdditionalDistributions.jl

@@ -67,6 +67,7 @@ module AdditionalDistributions
     include("continuous/CrystalBall.jl")
     include("continuous/Dagum.jl")
     include("continuous/Gompertz.jl")
+    include("continuous/IrwinHall.jl")
     include("continuous/Lomax.jl")
     include("continuous/Maxwell.jl")
     include("continuous/Nakagami.jl")
@@ -85,6 +86,7 @@ module AdditionalDistributions
         CrystalBall,
         Dagum,
         Gompertz,
+        IrwinHall,
         Lomax,
         Maxwell,
         Nakagami,

src/continuous/IrwinHall.jl

@@ -0,0 +1,158 @@
+
+"""
+    IrwinHall(n, a, b)
+
+``IrwinHall(n, a, b)`` distributed variable is a sum of ``n`` independent ``Uniform(a, b)`` variables. Its probability density function (PDF) is given by:
+
+```math
+f(x; n, a, b) = \\frac{1}{(b-a)(n-1)!}\\sum_{k=0}^{\\lfloor y\\rfloor}(-1)^k{n \\choose k} (y-k)^{n-1}, \\quad y = \\frac{x - a}{b - a}.
+```
+
+Related Bates distribution is ``IrwinHall(n, a/n, b/n)``.
+
+```julia
+IrwinHall(n)       # equivalent to IrwinHall(n, 0, 1)
+
+params(d)        # Get the parameters, i.e. (n, a, b)
+```
+
+External links:
+
+* [Irwin-Hall distribution on Wikipedia](https://en.wikipedia.org/wiki/Irwin%E2%80%93Hall_distribution)
+* [Bates distribution on Wikipedia](https://en.wikipedia.org/wiki/Bates_distribution)
+"""
+struct IrwinHall{S<:Integer, T<:Real} <: Distributions.ContinuousUnivariateDistribution
+    n::S
+    a::T
+    b::T
+    IrwinHall{S,T}(n, a, b) where {S<:Integer,T<:Real} = new{S,T}(n, a, b)
+end
+
+function IrwinHall(n::Integer, a::T, b::T; check_args::Bool=true) where {T<:Real}
+    @check_args IrwinHall (n, n > zero(n)) (b, b > a) 
+    return IrwinHall{typeof(n),typeof(a)}(n, a, b)
+end
+
+IrwinHall(n::Integer, a::Real, b::Real; check_args::Bool=true) = IrwinHall(n, promote(a, b)...; check_args=check_args)
+IrwinHall(n::Integer, a::Integer, b::Integer; check_args::Bool=true) = IrwinHall(n, float(a), float(b); check_args=check_args)
+
+IrwinHall(n::Integer) = IrwinHall(n, zero(n), one(n))
+
+@distr_support IrwinHall (d.n*d.a) (d.n*d.b)
+
+# parameters
+
+params(d::IrwinHall) = (d.n, d.a, d.b)
+
+shape(d::IrwinHall) = d.n
+location(d::IrwinHall) = d.n*d.a
+scale(d::IrwinHall) = d.b - d.a
+
+# moments, mode, median
+
+Statistics.mean(d::IrwinHall) = d.n*(d.a + d.b)/2
+Statistics.var(d::IrwinHall) = d.n*(d.b - d.a)^2 / 12
+StatsBase.skewness(d::IrwinHall) = zero(d.a)  
+StatsBase.kurtosis(d::IrwinHall) = -6 / (5*d.n)
+
+StatsBase.median(d::IrwinHall) = d.n*(d.a + d.b)/2
+StatsBase.mode(d::IrwinHall) = d.n*(d.a + d.b)/2
+
+# pdf, logpdf, cdf, quantile, cf, mgf, cgf
+
+function Distributions.pdf(d::IrwinHall, x::Real)
+    n, a, b = d.n, d.a, d.b
+    insupport(d, x) || return zero(a)
+
+    n == one(n) && return 1/(b - a) # uniform distribution
+
+    y = (x - n*a) / (b - a)
+    y = ifelse(y > n/2, n - y, y) # using symmetry we make calculation faster for large x
+    m = floor(Int, y)
+
+    c = one(y)/((b - a) * factorial(n - 1) )
+    S = c * y^(n - 1)
+
+    for k in 1:m
+        c *= -(n + 1 - k)/k
+        S += c * (y - k)^(n - 1)
+    end
+
+    return S
+end
+
+function Distributions.logpdf(d::IrwinHall, x::Real)
+    n, a, b = d.n, d.a, d.b
+    insupport(d, x) || return typemin(a)
+
+    y = (x - n*a) / (b - a)
+    y <= 1 && return (n-1)*log(y) - SpecialFunctions.logfactorial(n-1) - log(b - a) # left edge
+    y >= (n-1) && return (n-1)*log(n - y) - SpecialFunctions.logfactorial(n-1) - log(b - a) # right edge
+    return log(pdf(d, x)) # middle
+end
+
+function Distributions.cdf(d::IrwinHall, x::Real)
+    n, a, b = d.n, d.a, d.b
+    x <= minimum(d) && return zero(a)
+    x >= maximum(d) && return one(a)
+
+    n == one(n) && return (x - a)/(b - a) # uniform distribution
+
+    y = (x - n*a) / (b - a)
+    flag = y > n/2
+    y = ifelse(flag, n - y, y) # using symmetry we make calculation faster for large x
+    m = floor(Int, y)
+
+    c = one(y)/(n * factorial(n - 1) )
+    S = c * y^n
+
+    for k in 1:m
+        c *= -(n + 1 - k)/k
+        S += c * (y - k)^n
+    end
+
+    return flag ? 1-S : S
+end
+
+function Distributions.quantile(d::IrwinHall, p::Real)
+    n, a, b = d.n, d.a, d.b
+    A, B = minimum(d), maximum(d)
+    iszero(p) && return A
+    isone(p) && return B
+
+    p0 = 1/factorial(n)
+    p <= p0 && return A + (b-a) * (factorial(n)*p)^(1/n)
+    p >= 1-p0 && return B - (b-a) * (factorial(n)*(1-p))^(1/n)
+    return Roots.find_zero(x->cdf(d,x) - p, (A, B), Roots.Bisection())
+end
+
+Distributions.cf(d::IrwinHall, t::Real) = cf(Uniform(d.a, d.b), t) ^ d.n
+Distributions.mgf(d::IrwinHall, t::Real) = mgf(Uniform(d.a, d.b), t) ^ d.n
+Distributions.cgf(d::IrwinHall, t::Real) = d.n * cgf(Uniform(d.a, d.b), t)
+
+# affine transformations, convolution, conversion
+
+Base.:+(d::IrwinHall, x::Real) = IrwinHall(d.n, d.a + x/d.n, d.b + x/d.n)
+Base.:*(c::Real, d::IrwinHall) = IrwinHall(d.n, minmax(c*d.a, c*d.b)...)
+
+function Distributions.convolve(d1::IrwinHall, d2::IrwinHall)
+    d1.a ≈ d2.a || throw(ArgumentError("$(d1.a) !≈ $(d2.a): a parameters must be approximately equal"))
+    d1.b ≈ d2.b || throw(ArgumentError("$(d1.b) !≈ $(d2.b): b parameters must be approximately equal"))
+
+    return IrwinHall(d1.n + d2.n, d1.a, d1.b)
+end
+
+Distributions.convert(::Type{IrwinHall}, d::Uniform) = IrwinHall(1, d.a, d.b)
+
+function Distributions.convert(::Type{Uniform}, d::IrwinHall)
+    d.n != 1 && throw(DomainError("IrwinHall can be converted to Uniform only for n = 1"))
+    return Uniform(d.a, d.b)
+end
+
+function Distributions.convert(::Type{TriangularDist}, d::IrwinHall)
+    d.n != 2 && throw(DomainError("IrwinHall can be converted to TriangularDist only for n = 2"))
+    return TriangularDist(minimum(d), maximum(d), mean(d))
+end
+# sampling
+
+Distributions.rand(rng::Distributions.AbstractRNG, d::IrwinHall) = (d.b - d.a)*sum(rand(rng) for _ in 1:d.n) + d.n*d.a

test/continuous_test.jl

@@ -64,6 +64,12 @@ end
     M.check(Gompertz(0.1, 0.2))
 end
 
+@testitem "Generic – IrwinHall" tags=[:Generic, :Continuous, :IrwinHall] setup=[M] begin
+    M.check(IrwinHall(2))
+    M.check(IrwinHall(4, -1.0, 1.0))
+    M.check(IrwinHall(10, 2.0, 5.0))
+end
+
 @testitem "Generic – Lomax" tags=[:Generic, :Continuous, :Lomax] setup=[M] begin
     M.check(Lomax(0.5, 1.5))
     M.check(Lomax(2.0, 3.0))