Ekdist -- exponentiated Kumaraswamy

src/univariate/continuous/ekdist.jl

@@ -0,0 +1,188 @@
+"""
+    EKDist(α, β, γ)
+
+The *exponentiated Kumaraswamy (EK) distribution* with shape parameters `α > 0`, `β > 0`,
+and `γ > 0` has probability density function
+
+```math
+f(x; α, β, γ) = α β γ x^{α - 1} (x - x^α)^{β - 1} (1 - (1 - x^α)^β)^{γ - 1}, \\quad 0 < x < 1
+```
+
+When `γ == 1`, this reduces to the [Kumaraswamy distribution](@ref Kumaraswamy).
+
+References
+
+- Lemonte, A. J., Barreto-Souza, W., & Cordeiro, G. M. (2013). The exponentiated Kumaraswamy
+  distribution and its log-transform. Brazilian Journal of Probability and Statistics, 27(1),
+  31–53. http://www.jstor.org/stable/43601233
+"""
+struct EKDist{T<:Real} <: ContinuousUnivariateDistribution
+    α::T
+    β::T
+    γ::T
+    EKDist{T}(α, β, γ) where {T} = new{T}(α, β, γ)
+end
+
+function EKDist(α::T, β::T, γ::T; check_args::Bool=true) where {T<:Real}
+    @check_args EKDist (α, α > zero(α)) (β, β > zero(β)) (γ, γ > zero(γ))
+    return EKDist{T}(α, β, γ) 
+end
+
+EKDist(α::Real, β::Real, γ::Real; check_args::Bool=true) =
+    EKDist(float.(promote(α, β, γ))...; check_args=check_args)
+
+EKDist() = EKDist{Float64}(1.0, 1.0, 1.0)
+
+@distr_support EKDist 0 1
+
+#### Conversions
+convert(::Type{EKDist{T}}, α::S, β::S, γ::S, l::S, u::S) where {T <: Real, S <: Real} = 
+        EKDist(T(α), T(β), T(γ), T(l), T(u))
+
+convert(::Type{EKDist{T}}, d::EKDist) where {T<:Real} = 
+        EKDist{T}(T(d.α), T(d.β), T(d.γ), T(d.l), T(d.u))
+
+convert(::Type{EKDist{T}}, d::EKDist{T}) where {T<:Real} = d
+
+#### Parameters
+
+params(d::EKDist) = (d.α, d.β, d.γ, d.l, d.u)
+partype(::EKDist{T}) where {T} = T
+
+#### Evaluation
+
+function pdf(d::EKDist, x::Real)
+      xs = (x - d.l) / (d.u - d.l)
+      α = d.α
+      β = d.β
+      γ = d.γ
+    return ( α * β * γ * 
+             xs^(α - 1) * 
+             (1 - xs^α)^(β - 1) * 
+             (1 - (1 - xs^α)^β)^(γ - 1)
+           )
+end
+
+function logpdf(d::EKDist, x::Real)
+      xs =  (x - d.l) / (d.u - d.l) 
+      α = d.α
+      β = d.β
+      γ = d.γ
+    return  ( log(α) + log(β) + log(γ) + 
+              log(xs)*(α - 1) + 
+              log(1 - xs^α)*(β - 1) + 
+              log(1 - (1 - xs^α)^β)*(γ - 1)
+            ) 
+end
+
+function logpdf(d::EKDist, x::Real, edgebound::Real)
+    # Failes at exactly 1 or 0, so I guess jitter it by small amount
+    # this only matters to me in likelihood analysis for real data
+    # that might get defined as the max of the data or Fitting
+    # bound parameters as part of the process.
+    xc = clamp( x, d.l + edgebound, d.u - edgebound)
+    return logpdf(d, xc)  
+end
+
+# do we even need gradlogpdf?
+# grad implied gradient but this is definitely derivative wrt x
+function gradlogpdf(d::EKDist, x::Real)
+    outofsupport = (x < d.l) | (x > d.u)
+    xs = (x - d.l) / (d.u - d.l)
+    α = d.α
+    β = d.β
+    γ = d.γ
+
+    z = (α - 1) / xs +
+        -(α * xs^(α - 1)) / (1 - xs^α) * (β - 1) +
+        (α * β * xs^(α - 1)) * (1 - xs^α)^(β - 1) / (1 -(1 - xs^α)^β) * (γ - 1)
+
+    return outofsupport ? zero(z) : z
+end
+
+function gradlogpdf(d::EKDist, x::Real, edgebound::Real)
+    xc = clamp( x, d.l + edgebound, d.u - edgebound)
+    return gradlogpdf(d, xc)
+end
+
+
+
+function cdf(d::EKDist, x::Real)
+    α, β, γ = params(d)
+    y = (1 - (1 - clamp(x, 0, 1)^α)^β)^γ
+    return x < 0 ? zero(y) : (x > 1 ? one(y) : y)
+end
+
+function quantile(d::EKDist, q::Real)
+    α = d.α
+    β = d.β
+    γ = d.γ
+    xs = ( 1 - (1 - q^(1/γ) )^(1/β) )^(1/α) 
+    x = xs * (d.u - d.l) + d.l
+    return x
+end
+
+cquantile(d::EKDist, q::Real) = quantile(d, 1 - q)
+
+
+#### Sampling
+# `rand`: Uses fallback inversion sampling method
+
+## Fitting
+
+# `fit_mle` not yet implemented
+
+#### Statistics
+
+function ekdmoment(d::EKDist, k::Real)
+    α = d.α
+    β = d.β
+    γ = d.γ
+
+    # can't tell from paper if its infinite sum or stops at the floor of γ -1
+    # if infinite, we need to do convergence to get there and reflexive function on gamma(γ - i) when negative
+    if ( !isinteger(γ) | (γ < 0) ) return error("moment currently only implemented when γ is positive integer") end
+
+    # has infinite sums so probably need use Richardson.jl per googling
+
+     momentseries = [ β * gamma(γ+1) * (-1)^i * beta(1 + k/α, β * (i + 1) ) / (gamma(γ-i) * factorial(i)) for i in 0:Int(γ - 1)]
+
+     return sum(momentseries) 
+
+end
+
+mean(d::EKDist) = ekdmoment(d, 1)*(d.u - d.l) + d.l # need to scale it
+
+var(d::EKDist) = (ekdmoment(d, 2) - ekdmoment(d, 1)^2) * (d.u - d.l)^2 # no clue if this scaling makes sense...
+std(d::EKDist) = var(d::EKDist) ^.5
+
+# These aren't scaled for sure.
+skewness(d::EKDist) = if (d.l == 0) & (d.u == 1 )
+                        (ekdmoment(d, 3) - 3 * mean(d) * var(d) - mean(d)^3) / var(d)^(3/2) 
+                      else error("I didn't scale these yet, so it's not remotely right")
+                      end
+
+kurtosis(d::EKDist) = if (d.l == 0) & (d.u == 1 )
+                        a1 = ekdmoment(d, 1)
+                        a2 = ekdmoment(d, 2)
+                        a3 = ekdmoment(d, 3)
+                        a4 = ekdmoment(d, 4)
+
+                        return (a4 + a1 * (-4a3 + a1 * (6a2 - 3a1^2))) / (a2-a1^2)^2 - 3
+
+                    else return error("I didn't scale these yet, so it's not remotely right")
+                    end
+
+function mode(d::EKDist)
+    if (d.γ == 1) & (d.α >= 1) & (d.β >= 1) & !(d.α == d.β == 1)
+        return ( (d.α -1) / (d.α * d.β - 1) ) ^ (1 / d.α)
+
+    else (d.α > 1) & (d.β > 1) &  (d.γ > 1)
+        error("You'll have to approximate via pdf samples because I can't find general roots from derivatives like form: x^(α-1) * (1-x^α)^(β-1).
+                    \n Take caution because combinations of parameters can give you situations you  might not expect especially across the boundary of 1.0 for parameters. \n\n")
+
+    end
+end
+
+#= function entropy(d::EKDist) μ # In paper but not sure if is the expected form of entropy for package.
+end =#