diff --git a/docs/Project.toml b/docs/Project.toml index 540725c6d..2ef490a4d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" diff --git a/docs/src/univariate.md b/docs/src/univariate.md index 0b2c48c6e..a438f819f 100644 --- a/docs/src/univariate.md +++ b/docs/src/univariate.md @@ -139,6 +139,13 @@ Arcsine plotdensity((0.001, 0.999), Arcsine, (0, 1)) # hide ``` +```@docs +AsymmetricExponentialPower +``` +```@example plotdensity +plotdensity((-6, 8), AsymmetricExponentialPower, (0., 1., 2., 1., 0.7)) # hide +``` + ```@docs Beta ``` diff --git a/src/Distributions.jl b/src/Distributions.jl index ceb6063c7..b8b892927 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -64,6 +64,7 @@ export # distribution types Arcsine, + AsymmetricExponentialPower, Bernoulli, BernoulliLogit, Beta, diff --git a/src/univariate/continuous/asymmetricexponentialpower.jl b/src/univariate/continuous/asymmetricexponentialpower.jl new file mode 100644 index 000000000..8eb9dc817 --- /dev/null +++ b/src/univariate/continuous/asymmetricexponentialpower.jl @@ -0,0 +1,147 @@ +""" + AsymmetricExponentialPower(μ, σ, p₁, p₂, α) + +The *Asymmetric exponential power distribution*, with location `μ`, scale `σ`, left (`x < μ`) shape `p₁`, right (`x >= μ`) shape `p₂`, and skewness `α`, +has the probability density function [1] +```math +f(x; \\mu, \\sigma, p_1, p_2, \\alpha) = +\\begin{cases} +\\frac{\\alpha}{\\alpha^*} \\frac{1}{\\sigma} K_{EP}(p_1) \\exp \\left\\{ - \\frac{1}{2p_1}\\Big| \\frac{x-\\mu}{\\alpha^* \\sigma} \\Big|^{p_1} \\right\\}, & \\text{if } x \\leq \\mu \\\\ +\\frac{(1-\\alpha)}{(1-\\alpha^*)} \\frac{1}{\\sigma} K_{EP}(p_2) \\exp \\left\\{ - \\frac{1}{2p_2}\\Big| \\frac{x-\\mu}{(1-\\alpha^*) \\sigma} \\Big|^{p_2} \\right\\}, & \\text{if } x > \\mu +\\end{cases}, +``` +Where ``K_{EP}(p) = 1/(2p^{1/p}\\Gamma(1+1/p_1))`` and +```math +\\alpha^* = \\frac{\\alpha K_{EP}(p_1)}{\\alpha K_{EP}(p_1) + (1-\\alpha) K_{EP}(p_2)} +``` +The asymmetric exponential power distribution (AEPD) the skewed exponential power distribution as special case (``p_1 = p_2``) and thus the +Laplace, Normal, uniform, exponential power distribution, asymmetric Laplace and skew normal are also special cases. + +[1] Zhy, D. and V. Zinde-Walsh (2009). Properties and estimation of asymmetric exponential power distribution. _Journal of econometrics_, 148(1):86-96, 2009. + +```julia +AsymmetricExponentialPower() # AEPD with location 0, scale 1 left shape 2, right shape 2, and skewness 0.5 (the standard normal distribution) +AsymmetricExponentialPower(μ, σ, p₁, p₂ α) # AEPD with location μ, scale σ, left shape p₁, right shape p₂, and skewness α +AsymmetricExponentialPower(μ, σ, p₁, p₂) # AEPD with location μ, scale σ, left shape p₁, right shape p₂, and skewness 0.5 (the exponential power distribution) +AsymmetricExponentialPower(μ, σ) # AEPD with location μ, scale σ, left shape 2, right shape 2, and skewness 0.5 (the normal distribution) +AsymmetricExponentialPower(μ) # AEPD with location μ, scale 1, left shape 2, right shape 2, and skewness 0.5 (the normal distribution) + +params(d) # Get the parameters, i.e. (μ, σ, p₁, p₂, α) +shape(d) # Get the shape parameters, i.e. (p₁, p₂) +location(d) # Get the location parameter, i.e. μ +scale(d) # Get the scale parameter, i.e. σ +``` +""" +struct AsymmetricExponentialPower{T <: Real} <: ContinuousUnivariateDistribution + μ::T + σ::T + p₁::T + p₂::T + α::T + AsymmetricExponentialPower{T}(μ::T, σ::T, p₁::T, p₂::T, α::T) where {T} = new{T}(μ, σ, p₁, p₂, α) +end + + +function AsymmetricExponentialPower(µ::T, σ::T, p₁::T, p₂::T, α::T; check_args::Bool=true) where {T <: Real} + @check_args AsymmetricExponentialPower (σ, σ > zero(σ)) (p₁, p₁ > zero(p₁)) (p₂, p₂ > zero(p₂)) (α, zero(α) < α < one(α)) + return AsymmetricExponentialPower{T}(µ, σ, p₁, p₂, α) +end + +function AsymmetricExponentialPower(μ::Real, σ::Real, p₁::Real, p₂::Real, α::Real = 1//2; check_args::Bool=true) + return AsymmetricExponentialPower(promote(μ, σ, p₁, p₂, α)...; check_args=check_args) +end + +function AsymmetricExponentialPower(μ::Real, σ::Real; check_args::Bool=true) + return AsymmetricExponentialPower(promote(μ, σ, 2, 2, 1//2)...; check_args=check_args) +end + +AsymmetricExponentialPower(μ::Real=0) = AsymmetricExponentialPower(μ, 1, 2, 2, 1//2; check_args=false) + +@distr_support AsymmetricExponentialPower -Inf Inf + +### Conversions +function Base.convert(::Type{AsymmetricExponentialPower{T}}, d::AsymmetricExponentialPower) where {T<:Real} + AsymmetricExponentialPower{T}(T(d.μ), T(d.σ), T(d.p), T(d.α)) +end +Base.convert(::Type{AsymmetricExponentialPower{T}}, d::AsymmetricExponentialPower{T}) where {T<:Real} = d + +### Parameters +@inline partype(::AsymmetricExponentialPower{T}) where {T<:Real} = T +params(d::AsymmetricExponentialPower) = (d.μ, d.σ, d.p₁, d.p₂, d.α) +location(d::AsymmetricExponentialPower) = d.μ +shape(d::AsymmetricExponentialPower) = (d.p₁, d.p₂) +scale(d::AsymmetricExponentialPower) = d.σ + +# log of Equation (4) of [1] + +### Statistics +# Computes log K_{EP}(p), Zhy, D. and V. Zinde-Walsh (2009) +function logK(p::Real) + inv_p = inv(p) + return -(logtwo + loggamma(inv_p) + ((1 - p) * log(p))/p) +end + +# Equation (3) in Zhy, D. and V. Zinde-Walsh (2009) +function αstar(α::Real, p₁::Real, p₂::Real) + K1 = exp(logK(p₁)) + K2 = exp(logK(p₂)) + return α*K1 / (α*K1 + (1-α)*K2) +end + +# Computes Equation 4 in Zhy, D. and V. Zinde-Walsh (2009) +B(α::Real, p₁::Real, p₂::Real) = α*exp(logK(p₁)) + (1-α)*exp(logK(p₂)) + +#Calculates the kth central moment of the AEPD, Equation 14 in Zhy, D. and V. Zinde-Walsh (2009) +function m_k(d::AsymmetricExponentialPower, k::Integer) + _, σ, p₁, p₂, α = params(d) + inv_p1, inv_p2 = inv(p₁), inv(p₂) + H1 = k*log(p₁) + loggamma((1+k)*inv_p1) - (1+k)*loggamma(inv_p1) + H2 = k*log(p₂) + loggamma((1+k)*inv_p2) - (1+k)*loggamma(inv_p2) + return B(α, p₁, p₂)^(-k) * σ^k * ((-1)^k * α^(1+k)*exp(H1) + (1-α)^(1+k)*exp(H2)) +end + +mean(d::AsymmetricExponentialPower) = d.α == 1//2 ? float(d.μ) : m_k(d, 1) + d.μ +var(d::AsymmetricExponentialPower) = m_k(d, 2) - m_k(d, 1)^2 +skewness(d::AsymmetricExponentialPower) = d.α == 1//2 ? float(zero(partype(d))) : m_k(d, 3) / (std(d))^3 +kurtosis(d::AsymmetricExponentialPower) = m_k(d, 4)/var(d)^2 - 3 + +function logpdf(d::AsymmetricExponentialPower, x::Real) + μ, σ, p₁, p₂, α = params(d) + a, astar, inv_p, p = x < μ ? (α, αstar(α, p₁, p₂), inv(p₁), p₁) : (1 - α, 1-αstar(α, p₁, p₂), inv(p₂), p₂) + return -(log(astar) - log(a) + logtwo + log(σ) + loggamma(inv_p) + ((1 - p) * log(p) + (abs(μ - x) / (2 * σ * astar))^p) / p) +end + +function cdf(d::AsymmetricExponentialPower, x::Real) + μ, σ, p₁, p₂, α = params(d) + if x <= μ + inv_p = inv(p₁) + α * ccdf(Gamma(inv_p), inv_p * (abs((x-μ)/σ) / (2*αstar(α, p₁, p₂)))^p₁) + else + inv_p = inv(p₂) + α + (1-α) * cdf(Gamma(inv_p), inv_p * (abs((x-μ)/σ) / (2*(1-αstar(α, p₁, p₂))))^p₂) + end +end + +function quantile(d::AsymmetricExponentialPower, p::Real) + μ, σ, p₁, p₂, α = params(d) + if p <= α + inv_p = inv(p₁) + μ - 2*αstar(α, p₁, p₂)*σ * (p₁ * quantile(Gamma(inv_p), (α-p)/α))^inv_p + else + inv_p = inv(p₂) + μ + 2*(1-αstar(α, p₁, p₂))*σ * (p₂ * quantile(Gamma(inv_p), (p-α)/(1-α)))^inv_p + end +end + +function rand(rng::AbstractRNG, d::AsymmetricExponentialPower) + μ, σ, p₁, p₂, α = params(d) + if rand(rng) < α + inv_p = inv(p₁) + z = 2*σ * (p₁ * rand(rng,Gamma(inv_p, 1)))^inv_p + return μ - αstar(α,p₁,p₂) * z + else + inv_p = inv(p₂) + z = 2*σ * (p₂ * rand(rng,Gamma(inv_p, 1)))^inv_p + return μ + (1-αstar(α,p₁,p₂)) * z + end +end \ No newline at end of file diff --git a/src/univariates.jl b/src/univariates.jl index 0e2ae05e3..754db31cd 100644 --- a/src/univariates.jl +++ b/src/univariates.jl @@ -669,6 +669,7 @@ const discrete_distributions = [ const continuous_distributions = [ "arcsine", + "asymmetricexponentialpower", "beta", "betaprime", "biweight", diff --git a/test/runtests.jl b/test/runtests.jl index ce3f16b79..1a24bc518 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,6 +10,8 @@ using LinearAlgebra import JSON import ForwardDiff + + const tests = [ "univariate/continuous/loguniform", "univariate/continuous/arcsine", @@ -87,6 +89,7 @@ const tests = [ "density_interface", "reshaped", "univariate/continuous/skewedexponentialpower", + "univariate/continuous/asymmetricexponentialpower", "univariate/discrete/discreteuniform", "univariate/continuous/tdist", "univariate/orderstatistic", @@ -154,7 +157,7 @@ const tests = [ # "univariate/discrete/geometric", # "univariate/discrete/hypergeometric", # "univariate/discrete/noncentralhypergeometric", - # "univariate/discrete/poisson", + # "univariate/discrete/poisson",a # "univariate/discrete/skellam", ### file is present but was not included in list diff --git a/test/univariate/continuous/asymmetricexponentialpower.jl b/test/univariate/continuous/asymmetricexponentialpower.jl new file mode 100644 index 000000000..53a444540 --- /dev/null +++ b/test/univariate/continuous/asymmetricexponentialpower.jl @@ -0,0 +1,90 @@ +using Distributions +using SpecialFunctions + +using Test + +@testset "AsymmetricExponentialPower" begin + @testset "p₁ = p₂" begin + @testset "α = 0.5" begin + @test_throws DomainError AsymmetricExponentialPower(0, 0, 0, 0, 0) + d1 = AsymmetricExponentialPower(0, 1, 1, 1, 0.5f0) + @test @inferred partype(d1) == Float32 + d2 = AsymmetricExponentialPower(0, 1, 1, 1, 0.5) + @test @inferred partype(d2) == Float64 + @test @inferred params(d2) == (0., 1., 1., 1., 0.5) + @test @inferred cdf(d2, Inf) == 1 + @test @inferred cdf(d2, -Inf) == 0 + @test @inferred quantile(d2, 1) == Inf + @test @inferred quantile(d2, 0) == -Inf + test_distr(d2, 10^6) + + # Comparison to SEPD + d = AsymmetricExponentialPower(0, 1, 1, 1, 0.5) + dl = SkewedExponentialPower(0, 1, 1, 0.5) + @test @inferred mean(d) ≈ mean(dl) + @test @inferred var(d) ≈ var(dl) + @test @inferred skewness(d) ≈ skewness(dl) + @test @inferred kurtosis(d) ≈ kurtosis(dl) + @test @inferred pdf(d, 0.5) ≈ pdf(dl, 0.5) + @test @inferred cdf(d, 0.5) ≈ cdf(dl, 0.5) + @test @inferred quantile(d, 0.5) ≈ quantile(dl, 0.5) + + # Comparison to laplace + d = AsymmetricExponentialPower(0, 1, 1, 1, 0.5) + dl = Laplace(0, 1) + @test @inferred mean(d) ≈ mean(dl) + @test @inferred var(d) ≈ var(dl) + @test @inferred skewness(d) ≈ skewness(dl) + @test @inferred kurtosis(d) ≈ kurtosis(dl) + @test @inferred pdf(d, 0.5) ≈ pdf(dl, 0.5) + @test @inferred cdf(d, 0.5) ≈ cdf(dl, 0.5) + @test @inferred quantile(d, 0.5) ≈ quantile(dl, 0.5) + + # comparison to exponential power distribution (PGeneralizedGaussian), + # where the variance is reparametrized as σₚ = p^(1/p)σ to ensure equal pdfs + p = 1.2 + d = AsymmetricExponentialPower(0, 1, p, p, 0.5) + de = PGeneralizedGaussian(0, p^(1/p), p) + @test @inferred mean(d) ≈ mean(de) + @test @inferred var(d) ≈ var(de) + @test @inferred skewness(d) ≈ skewness(de) + @test @inferred kurtosis(d) ≈ kurtosis(de) + @test @inferred pdf(d, 0.5) ≈ pdf(de, 0.5) + @test @inferred cdf(d, 0.5) ≈ cdf(de, 0.5) + test_distr(d, 10^6) + end + + @testset "α != 0.5" begin + # relationship between aepd(μ, σ, p, p, α) and + # aepd(μ, σ, p, p, 1-α) + d1 = AsymmetricExponentialPower(0, 1, 0.1, 0.1, 0.7) + d2 = AsymmetricExponentialPower(0, 1, 0.1, 0.1, 0.3) + @inferred -mean(d1) ≈ mean(d2) + @inferred var(d1) ≈ var(d2) + @inferred -skewness(d1) ≈ skewness(d2) + @inferred kurtosis(d1) ≈ kurtosis(d2) + + α, p = rand(2) + d = SkewedExponentialPower(0, 1, p, α) + # moments of the standard SEPD, Equation 18 in Zhy, D. and V. Zinde-Walsh (2009) + moments = [(2*p^(1/p))^k * ((-1)^k*α^(1+k) + (1-α)^(1+k)) * gamma((1+k)/p)/gamma(1/p) for k ∈ 1:4] + + @inferred var(d) ≈ moments[2] - moments[1]^2 + @inferred skewness(d) ≈ moments[3] / (√(moments[2] - moments[1]^2))^3 + @inferred kurtosis(d) ≈ (moments[4] / ((moments[2] - moments[1]^2))^2 - 3) + test_distr(d, 10^6) + end + end + @testset "p₁ ≠ p₂" begin + d2 = AsymmetricExponentialPower(0, 1., 1., 1., 0.7) + test_distr(d2, 10^6) + + # relationship between aepd(0, σ, p₁, p₂, α) and aepd(0, σ, p₂, p₁, 1-α) + d1 = AsymmetricExponentialPower(0, 1., 0.4, 1.2, 0.7) + d2 = AsymmetricExponentialPower(0, 1., 1.2, 0.4, 0.3) + @inferred -mean(d1) ≈ mean(d2) + @inferred var(d1) ≈ var(d2) + @inferred -skewness(d1) ≈ skewness(d2) + @inferred kurtosis(d1) ≈ kurtosis(d2) + end +end \ No newline at end of file