diff --git a/docs/src/univariate.md b/docs/src/univariate.md index 3c2e6d4dc..f8ee6f164 100644 --- a/docs/src/univariate.md +++ b/docs/src/univariate.md @@ -244,6 +244,20 @@ GeneralizedExtremeValue plotdensity((0, 30), GeneralizedExtremeValue, (0, 1, 1)) # hide ``` +```@docs +GeneralizedHyperbolic +``` +```@example plotdensity +plotdensity((-8, 8), GeneralizedHyperbolic, (3, 2, 1, -1, 1/2)) # hide +``` + +```@docs +GeneralizedInverseGaussian +``` +```@example plotdensity +plotdensity((0, 20), GeneralizedInverseGaussian, (3, 8, -1/2)) # hide +``` + ```@docs GeneralizedPareto ``` diff --git a/src/Distributions.jl b/src/Distributions.jl index fcf699434..39574386a 100644 --- a/src/Distributions.jl +++ b/src/Distributions.jl @@ -100,6 +100,8 @@ export FullNormalCanon, Gamma, DiscreteNonParametric, + GeneralizedHyperbolic, + GeneralizedInverseGaussian, GeneralizedPareto, GeneralizedExtremeValue, Geometric, diff --git a/src/univariate/continuous/generalizedhyperbolic.jl b/src/univariate/continuous/generalizedhyperbolic.jl new file mode 100644 index 000000000..9dbbd71e0 --- /dev/null +++ b/src/univariate/continuous/generalizedhyperbolic.jl @@ -0,0 +1,328 @@ +@doc raw""" + GeneralizedHyperbolic(α, β, δ, μ=0, λ=1) + +The *generalized hyperbolic (GH) distribution* with traditional parameters: + +- $\alpha>0$ (shape); +- $-\alpha<\beta<\alpha$ (skewness); +- $\delta>0$ ("scale", but not really, because it appears as an argument to the modified Bessel function of the 2nd kind in the normalizing constant); +- $\mu\in\mathbb R$ (location); +- $\lambda\in\mathbb R$ is a shape parameter, where $\lambda\neq 1$ makes the distribution "generalized" + +has probability density function: + +```math +\frac{ + (\gamma/\delta)^{\lambda} +}{ + \sqrt{2\pi} K_{\lambda}(\delta \gamma) +} +e^{\beta (x-\mu)} +\frac{ + K_{\lambda-1/2}\left(\alpha\sqrt{\delta^2 + (x-\mu)^2}\right) +}{ + \left(\alpha^{-1} \sqrt{\delta^2 + (x-\mu)^2}\right)^{1/2 - \lambda} +}, \quad\gamma=\sqrt{\alpha^2 - \beta^2} +``` + +These parameters are actually stored in `struct GeneralizedHyperbolic{T<:Real}`. + +External links: + +* [Generalized hyperbolic distribution on Wikipedia](https://en.wikipedia.org/wiki/Generalised_hyperbolic_distribution). +""" +struct GeneralizedHyperbolic{T<:Real} <: ContinuousUnivariateDistribution + α::T + β::T + δ::T + μ::T + λ::T + function GeneralizedHyperbolic(α::T, β::T, δ::T, μ::T=zero(T), λ::T=one(T); check_args::Bool=true) where T<:Real + check_args && @check_args GeneralizedHyperbolic (α, α > zero(α)) (δ, δ > zero(δ)) (β, -α < β < α) + new{T}(α, β, δ, μ, λ) + end +end + +GeneralizedHyperbolic(α::Real, β::Real, δ::Real, μ::Real=0, λ::Real=1; check_args::Bool=true) = + GeneralizedHyperbolic(promote(α, β, δ, μ, λ)...; check_args) + +@doc raw""" + GeneralizedHyperbolic(Val(:locscale), z, p=0, μ=0, σ=1, λ=1) + +Location-scale parameterization [1] of the generalized hyperbolic distribution with parameters + +- $z>0$ (shape); +- $p\in\mathbb R$ measures skewness ($p=0$ results in a symmetric distribution); +- $\mu\in\mathbb R$ and $\sigma>0$ are location and scale; +- $\lambda\in\mathbb R$ is a shape parameter, where $\lambda\neq 1$ makes the distribution "generalized" + +has probability density function: + +```math +\frac{\sqrt z}{ + \sqrt{2\pi} K_{\lambda}(z) +} +e^{p z \cdot\varepsilon} +\sqrt{ + \left(\frac{1+\varepsilon^2}{1+p^2}\right)^{\lambda - 1/2} +} +K_{\lambda-1/2}\left[ + z \sqrt{(1+p^2)(1+\varepsilon^2)} +\right] +``` + +These parameters are _not_ stored in `struct GeneralizedHyperbolic`. +Use `params(d, Val(:locscale))`, where `d` is an instance of `GeneralizedHyperbolic`, to retrieve them. + +Advantages of this parameterization: + +- It's truly location-scale, whereas $\delta$ in the traditional parameterization isn't a true scale parameter. +- All parameters are either positive or unconstrained. The traditional parameterization has the complicated linear constraint $-\alpha<\beta<\alpha$. + +References: + +1. Puig, Pedro, and Michael A. Stephens. “Goodness-of-Fit Tests for the Hyperbolic Distribution.” The Canadian Journal of Statistics / La Revue Canadienne de Statistique 29, no. 2 (2001): 309–20. https://doi.org/10.2307/3316079. +""" +GeneralizedHyperbolic(::Val{:locscale}, z::Real, p::Real=0, μ::Real=0, σ::Real=1, λ::Real=1; check_args::Bool=true) = + GeneralizedHyperbolic(z * sqrt(1 + p^2)/σ, z * p / σ, σ, μ, λ; check_args) + +params(d::GeneralizedHyperbolic) = (d.α, d.β, d.δ, d.μ, d.λ) +params(d::GeneralizedHyperbolic, ::Val{:locscale}) = begin + α, β, δ, μ, λ = params(d) + γ = sqrt(α^2 - β^2) + + (; z=δ * γ, p=β / γ, μ, σ=δ, λ) +end +partype(::GeneralizedHyperbolic{T}) where T = T + +minimum(::GeneralizedHyperbolic) = -Inf +maximum(::GeneralizedHyperbolic) = Inf +insupport(::GeneralizedHyperbolic, x::Real) = true + +"Fit quadratic `y = ax^2 + bx + c` through 3 points (x, y), return coefficients `(a, b)`." +function _fit_quadratic(x1, y1, x2, y2, x3, y3) + @assert x1 <= x2 <= x3 + denom = (x1-x2) * (x1-x3) * (x2-x3) # < 0 + a = ( + x3 * (y2-y1) + x2 * (y1-y3) + x1 * (y3-y2) + ) / denom + b = ( + x3^2 * (y1-y2) + x1^2 * (y2-y3) + x2^2 * (y3-y1) + ) / denom + (a, b) +end + +""" + mode(::GeneralizedHyperbolic) + +- Exact formulae are used for λ=1 and λ=2. +- For other values of λ quadratic fit search is used. The initial bracket `lo < mode < hi` is computed based on +inequalities from [1, eq. 2.27]. + +## References + +1. Robert E. Gaunt, Milan Merkle, "On bounds for the mode and median of the generalized hyperbolic and related distributions", Journal of Mathematical Analysis and Applications, Volume 493, Issue 1, 2021, 124508, ISSN 0022-247X, https://doi.org/10.1016/j.jmaa.2020.124508. +""" +function mode(d::GeneralizedHyperbolic) + α, β, δ, μ, λ = params(d) + γ = sqrt(α^2 - β^2) + + if λ ≈ 1 + μ + β * δ / γ # Wolfram + elseif λ ≈ 2 + μ + β / α / γ^2 * (α + sqrt(β^2 + (α * δ * γ)^2)) # Wolfram + else + lo, hi = let # Bounds for the bracketing interval + EX = mean(d) + # eq. 2.27 from the paper: EX - upper < mode < EX - lower for β>0. + # When β<0, the inequality is reversed. + lower = β/γ^2 * ( + 1/2 + sqrt(λ^2 + δ^2 * γ^2) - sqrt((λ - 1/2)^2 + δ^2 * γ^2) + ) + upper = β/γ^2 * ( + 5/2 + sqrt((λ + 1)^2 + δ^2 * γ^2) - sqrt((λ - 3/2)^2 + δ^2 * γ^2) + ) + if β < 0 + upper, lower = lower, upper + end + EX - upper, EX - lower + end + + # Minimize negative log-PDF using quadratic fit search + x1, x2, x3 = lo, (lo + hi)/2, hi # invariant: x1 < x2 < x3 + xopt = x2 + y1, y2, y3 = -logpdf(d, x1), -logpdf(d, x2), -logpdf(d, x3) + a, b = _fit_quadratic(x1, y1, x2, y2, x3, y3) + if a < 0 + @warn "Quadratic points up instead of down: a=$a." + return xopt + end + (!isfinite(a) || !isfinite(b)) && return xopt + niter = 0 + while x3 - x1 > 1e-6 + niter += 1 + xopt = -b / (2a) + yopt = -logpdf(d, xopt) + + if xopt < x2 + if yopt > y2 + x1, y1 = xopt, yopt + else + x2, x3 = xopt, x2 + y2, y3 = yopt, y2 + end + else # xopt > x2 + if yopt > y2 + x3, y3 = xopt, yopt + else + x1, x2 = x2, xopt + y1, y2 = y2, yopt + end + end + + a, b = _fit_quadratic(x1, y1, x2, y2, x3, y3) + if !isfinite(a) || !isfinite(b) + (abs(x2 - x1) < 1e-6 || abs(x3 - x2) < 1e-6) && break + @warn "Failed to build quadratic" (; niter, x1, y1, x2, y2, x3, y3) + break + end + end + xopt + end +end + +mean(d::GeneralizedHyperbolic) = begin + α, β, δ, μ, λ = params(d) + γ = sqrt(α^2 - β^2) + μ + β * δ / γ * besselk(1 + λ, δ * γ) / besselk(λ, δ * γ) +end + +var(d::GeneralizedHyperbolic) = begin + α, β, δ, μ, λ = params(d) + γ = sqrt(α^2 - β^2) + + t0 = besselk(0 + λ, δ * γ) + t1 = besselk(1 + λ, δ * γ) + t2 = besselk(2 + λ, δ * γ) + δ / γ * t1/t0 - (β * δ / γ * t1/t0)^2 + (β * δ / γ)^2 * t2/t0 +end + +skewness(d::GeneralizedHyperbolic) = begin + α, β, δ, μ, λ = params(d) + γ = sqrt(α^2 - β^2) + + t0 = besselk(0 + λ, δ * γ) + t1 = besselk(1 + λ, δ * γ) + t2 = besselk(2 + λ, δ * γ) + t3 = besselk(3 + λ, δ * γ) + ( + -3β * (δ / γ * t1/t0)^2 + 2 * (β * δ / γ * t1/t0)^3 + 3β * (δ / γ)^2 * t2/t0 + - 3 * (β * δ / γ)^3 * t1*t2/t0^2 + (β * δ / γ)^3 * t3/t0 + ) / sqrt( + δ / γ * t1/t0 - (β * δ / γ * t1/t0)^2 + (β * δ / γ)^2 * t2/t0 + )^3 +end + +kurtosis(d::GeneralizedHyperbolic) = begin + α, β, δ, μ, λ = params(d) + γ = sqrt(α^2 - β^2) + + t0 = besselk(0 + λ, δ * γ) + t1 = besselk(1 + λ, δ * γ) + t2 = besselk(2 + λ, δ * γ) + t3 = besselk(3 + λ, δ * γ) + t4 = besselk(4 + λ, δ * γ) + ( + 3 * γ^2 * t0^3 * t2 + 6 * β^2 * γ * δ * t0 * (t1^3 - 2t0 * t1 * t2 + t0^2 * t3) + + β^4 * δ^2 * (-3 * t1^4 + 6t0 * t1^2 * t2 - 4 * t0^2 * t1 * t3 + t0^3 * t4) + ) / ( + γ * t0 * t1 + β^2 * δ * (-t1^2 + t0 * t2) + )^2 - 3 # EXCESS kurtosis +end + +logpdf(d::GeneralizedHyperbolic, x::Real) = begin + α, β, δ, μ, λ = params(d) + γ = sqrt(α^2 - β^2) + + function logbesselk(v::Real, x::Real, K::Integer=5) + if x > 600 + # Asymptotic expansion, works for massive values of `x` on the order of 10^3, 10^4 and higher. + # Important, because otherwise PDF becomes exactly zero `exp(-Inf)==0` way too early. + μ = 4v^2 + term = one(x) + s = one(x) + for k in 1:K + term *= (μ - (2k-1)^2) / (k * 8x) + s += term + end + (log(π) - log(2x))/2 - x + log(abs(s)) + else + log(besselk(v, x)) # Returns `-Inf` for x>600 + end + end + + ( + -0.5log(2π) - logbesselk(λ, γ * δ) + λ * (log(γ) - log(δ)) + + β * (x - μ) + + (λ - 1/2) * (0.5log(δ^2 + (x - μ)^2) - log(α)) + + logbesselk(λ - 1/2, α * sqrt(δ^2 + (x - μ)^2)) + ) +end + +cdf(d::GeneralizedHyperbolic, x::Real) = + if isinf(x) + (x < 0) ? zero(x) : one(x) + elseif isnan(x) + typeof(x)(NaN) + else + quadgk(z -> pdf(d, z), -Inf, x, maxevals=10^4)[1] + end + +@quantile_newton GeneralizedHyperbolic + +mgf(d::GeneralizedHyperbolic, t::Number) = begin + α, β, δ, μ, λ = params(d) + γ = sqrt(α^2 - β^2) + + g = sqrt(α^2 - (t + β)^2) + exp(t * μ) / g^λ * sqrt((α - β) * (α + β))^λ * besselk(λ, g * δ) / besselk(λ, γ * δ) +end + +cf(d::GeneralizedHyperbolic, t::Number) = mgf(d, 1im * t) + +@doc raw""" + rand(::AbstractRNG, ::GeneralizedHyperbolic) + +Sample from `GeneralizedHyperbolic(α, β, δ, μ, λ)` using its mixture representation: + +```math +\begin{aligned} +\gamma &= \sqrt{\alpha^2 - \beta^2}\\ +V &\sim \mathrm{GeneralizedInverseGaussian}(\delta / \gamma, \delta^2, \lambda)\\ +\xi &= \mu + \beta V + \sqrt{V} \varepsilon, \quad\varepsilon \sim \mathcal N(0,1) +\end{aligned} +``` + +Then ξ is distributed as `GeneralizedHyperbolic(α, β, δ, μ, λ)`. + +Verified in Wolfram Mathematica: + +``` +In:= TransformedDistribution[\[Mu] + \[Beta]*V + + Sqrt[V] \[Epsilon], {\[Epsilon] \[Distributed] NormalDistribution[], + V \[Distributed] + InverseGaussianDistribution[\[Delta]/ + Sqrt[\[Alpha]^2 - \[Beta]^2], \[Delta]^2, \[Lambda]]}] + +Out= HyperbolicDistribution[\[Lambda], \[Alpha], \[Beta], \[Delta], \[Mu]] +``` + +Note that here λ is the first parameter, while in this implementation it's the _last_ one. +""" +rand(rng::AbstractRNG, d::GeneralizedHyperbolic) = begin + α, β, δ, μ, λ = params(d) + γ = sqrt(α^2 - β^2) + + V = rand(rng, GeneralizedInverseGaussian(δ/γ, δ^2, λ)) + μ + β * V + sqrt(V) * randn(rng) +end diff --git a/src/univariate/continuous/generalizedinversegaussian.jl b/src/univariate/continuous/generalizedinversegaussian.jl new file mode 100644 index 000000000..f5f6c9e9b --- /dev/null +++ b/src/univariate/continuous/generalizedinversegaussian.jl @@ -0,0 +1,242 @@ +import SpecialFunctions: besselk + +@doc raw""" + GeneralizedInverseGaussian(μ, λ, θ) + +The *generalized inverse Gaussian distribution* with parameters `\mu>0`, `\lambda>0` and real `\theta` has probability density function: + +```math +f(x; \mu, \lambda, \theta) = +\frac{1}{ + 2\mu^{\theta} K_{\theta}(\lambda/\mu) +} x^{\theta-1} \exp\left( + -\frac{\lambda}{2} \left(\frac{1}{x} + \frac{x}{\mu^2}\right) +\right) +, \quad x > 0 +``` + +External links: + +* [Generalized Inverse Gaussian distribution on Wikipedia](https://en.wikipedia.org/wiki/Generalized_inverse_Gaussian_distribution). +""" +struct GeneralizedInverseGaussian{T<:Real} <: ContinuousUnivariateDistribution + μ::T + λ::T + θ::T + function GeneralizedInverseGaussian(μ::T, λ::T, θ::T; check_args::Bool=true) where T<:Real + check_args && @check_args GeneralizedInverseGaussian (μ, μ > zero(μ)) (λ, λ > zero(λ)) + new{T}(μ, λ, θ) + end +end + +GeneralizedInverseGaussian(μ::Real, λ::Real, θ::Real=-1/2; check_args::Bool=true) = + GeneralizedInverseGaussian(promote(μ, λ, θ)...; check_args) + +params(d::GeneralizedInverseGaussian) = (d.μ, d.λ, d.θ) +partype(::GeneralizedInverseGaussian{T}) where T = T + +minimum(::GeneralizedInverseGaussian) = 0.0 +maximum(::GeneralizedInverseGaussian) = Inf +insupport(::GeneralizedInverseGaussian, x::Real) = x >= 0 + +mode(d::GeneralizedInverseGaussian) = begin + μ, λ, θ = params(d) + tmp = μ * (θ - 1) + μ / λ * ( + tmp + sqrt(λ^2 + tmp^2) + ) +end + +mean(d::GeneralizedInverseGaussian) = begin + μ, λ, θ = params(d) + b0 = besselk(0 + θ, λ / μ) + b1 = besselk(1 + θ, λ / μ) + μ * b1 / b0 +end + +var(d::GeneralizedInverseGaussian) = begin + μ, λ, θ = params(d) + b0 = besselk(0 + θ, λ / μ) + b1 = besselk(1 + θ, λ / μ) + b2 = besselk(2 + θ, λ / μ) + + μ^2 * (b0 * b2 - b1^2) / b0^2 +end + +# Source: Wolfram +skewness(d::GeneralizedInverseGaussian) = begin + μ, λ, θ = params(d) + b0 = besselk(0 + θ, λ / μ) + b1 = besselk(1 + θ, λ / μ) + b2 = besselk(2 + θ, λ / μ) + b3 = besselk(3 + θ, λ / μ) + ( + 2 * b1^3 - 3b0 * b1 * b2 + b0^2 * b3 + ) / sqrt( + b0 * b2 - b1^2 + )^3 +end + +# Source: Wolfram +kurtosis(d::GeneralizedInverseGaussian) = begin + μ, λ, θ = params(d) + t0 = besselk(0 + θ, λ/μ) + t1 = besselk(1 + θ, λ/μ) + t2 = besselk(2 + θ, λ/μ) + t3 = besselk(3 + θ, λ/μ) + t4 = besselk(4 + θ, λ/μ) + ( + -3 * t1^4 + 6t0 * t1^2 * t2 - 4 * t0^2 * t1 * t3 + t0^3 * t4 + ) / ( + t1^2 - t0 * t2 + )^2 - 3 # EXCESS kurtosis! +end + +logpdf(d::GeneralizedInverseGaussian, x::Real) = begin + μ, λ, θ = params(d) + if x >= 0 + -log(2) - θ * log(μ) - log(besselk(θ, λ / μ)) - λ/2 * (1/x + x/μ^2) + (θ - 1) * log(x) + else + -Inf + end +end + +cdf(d::GeneralizedInverseGaussian, x::Real) = + if isinf(x) + (x < 0) ? zero(x) : one(x) + elseif isnan(x) + typeof(x)(NaN) + else + quadgk(z -> pdf(d, z), 0, x, maxevals=1000)[1] + end + +@quantile_newton GeneralizedInverseGaussian + +mgf(d::GeneralizedInverseGaussian, t::Number) = begin + μ, λ, θ = params(d) + tmp = 1 - 2t * μ^2 / λ + tmp^(-θ/2) * besselk(θ, λ/μ * sqrt(tmp)) / besselk(θ, λ/μ) +end + +cf(d::GeneralizedInverseGaussian, t::Number) = mgf(d, 1im * t) + +""" + rand(rng::AbstractRNG, d::GeneralizedInverseGaussian) + +Sample from the generalized inverse Gaussian distribution based on [1], end of Section 6. + +### References + +1. Devroye, Luc. 2014. “Random Variate Generation for the Generalized Inverse Gaussian Distribution.” Statistics and Computing 24 (2): 239–46. https://doi.org/10.1007/s11222-012-9367-z. +""" +rand(rng::AbstractRNG, d::GeneralizedInverseGaussian) = begin + # If X ~ GIG(μ, λ, θ), then X/μ ~ _GIG(λ/μ, θ), + # so mu * _GIG(λ/μ, θ) == GIG(μ, λ, θ) + μ, λ, θ = params(d) + μ * rand(rng, _GIG(λ/μ, θ)) +end + +# ===== Private two-parameter version ===== +""" + _GIG(λ, ω) + +Two-parameter generalized inverse Gaussian distribution, only used for sampling. + +If `X ~ GeneralizedInverseGaussian(μ, λ, θ)`, then `Y = X / μ` follows `_GIG(λ/μ, θ)`. +NOTE: the paper says (Section 1) that the second parameter of `_GIG` should be `ω = 2sqrt(b/a)`, but computations in Wolfram Mathematica show otherwise. + +### References + +1. Devroye, Luc. 2014. “Random Variate Generation for the Generalized Inverse Gaussian Distribution.” Statistics and Computing 24 (2): 239–46. https://doi.org/10.1007/s11222-012-9367-z. +""" +struct _GIG{T1<:Real, T2<:Real} <: ContinuousUnivariateDistribution + ω::T1 + λ::T2 + function _GIG(ω::T1, λ::T2) where {T1<:Real, T2<:Real} + @assert ω >= 0 + new{T1, T2}(ω, λ) + end +end + +logpdf(d::_GIG, x::Real) = + if x >= 0 + -log(2 * besselk(d.λ, d.ω)) + (d.λ - 1) * log(x) - d.ω/2 * (x + 1/x) + else + -Inf + end + +""" + rand(rng::AbstractRNG, d::_GIG) + +Sampling from the _2-parameter_ generalized inverse Gaussian distribution based on [1], end of Section 6. + +### References + +1. Devroye, Luc. 2014. “Random Variate Generation for the Generalized Inverse Gaussian Distribution.” Statistics and Computing 24 (2): 239–46. https://doi.org/10.1007/s11222-012-9367-z. +""" +function rand(rng::AbstractRNG, d::_GIG) + λ, ω = d.λ, d.ω + negative_lmb = λ < 0 + λ = abs(λ) + + α = sqrt(ω^2 + λ^2) - λ + ψ(x) = -α * (cosh(x) - 1) - λ * (exp(x) - x - 1) + ψprime(x) = -α * sinh(x) - λ * (exp(x) - 1) + + tmp = -ψ(1) + t = if 0.5 <= tmp <= 2 + 1.0 + elseif tmp > 2 + sqrt(2 / (α + λ)) + else + log(4 / (α + 2λ)) + end + + tmp = -ψ(-1) + s = if 0.5 <= tmp <= 2 + 1.0 + elseif tmp > 2 + sqrt(4 / (α * cosh(1) + λ)) + else + min(1/λ, log(1 + 1/α + sqrt(1 / α^2 + 2/α))) + end + + eta, zeta, theta, xi = -ψ(t), -ψprime(t), -ψ(-s), ψprime(-s) + p, r = 1/xi, 1/zeta + + t_ = t - r * eta + s_ = s - p * theta + q = t_ + s_ + + chi(x) = if -s_ <= x <= t_ + 1.0 + elseif x < -s_ + exp(-theta + xi * (x + s)) + else # x > t_ + exp(-eta - zeta * (x - t)) + end + + # Generation + U, V, W = rand(rng), rand(rng), rand(rng) + X = if U < q / (p + q + r) + -s_ + q * V + elseif U < (q + r) / (p + q + r) + t_ - r * log(V) + else + -s_ + p * log(V) + end + while W * chi(X) > exp(ψ(X)) + U, V, W = rand(rng), rand(rng), rand(rng) + X = if U < q / (p + q + r) + -s_ + q * V + elseif U < (q + r) / (p + q + r) + t_ - r * log(V) + else + -s_ + p * log(V) + end + end + + tmp = λ/ω + logGIG = log(tmp + sqrt(1 + tmp^2)) + X + exp(negative_lmb ? -logGIG : logGIG) +end diff --git a/src/univariates.jl b/src/univariates.jl index b60e5a294..d2e0f94cc 100644 --- a/src/univariates.jl +++ b/src/univariates.jl @@ -693,8 +693,10 @@ const continuous_distributions = [ "frechet", "gamma", "erlang", "pgeneralizedgaussian", # GeneralizedGaussian depends on Gamma + "generalizedhyperbolic", "generalizedpareto", "generalizedextremevalue", + "generalizedinversegaussian", "gumbel", "inversegamma", "inversegaussian", diff --git a/test/runtests.jl b/test/runtests.jl index f63a3d5e8..8ef2fc8a0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,6 +15,7 @@ const tests = [ "univariate/continuous/loguniform", "univariate/continuous/arcsine", "univariate/discrete/dirac", + "univariate/continuous/generalizedhyperbolic", "truncate", "truncated/normal", "truncated/exponential", @@ -100,6 +101,7 @@ const tests = [ "univariate/continuous/triangular", "statsapi", "univariate/continuous/inversegaussian", + "univariate/continuous/generalizedinversegaussian", ### missing files compared to /src: # "common", diff --git a/test/univariate/continuous/generalizedhyperbolic.jl b/test/univariate/continuous/generalizedhyperbolic.jl new file mode 100644 index 000000000..83bfb4fe4 --- /dev/null +++ b/test/univariate/continuous/generalizedhyperbolic.jl @@ -0,0 +1,175 @@ +@testset "Generalized hyperbolic" begin + Hyp(z, p=0, μ=0, σ=1, λ=1) = GeneralizedHyperbolic(Val(:locscale), z, p, μ, σ, λ) + # Empirical characteristic function + cf_empirical(samples::AbstractVector{<:Real}, t::Real) = mean(x->exp(1im * t * x), samples) + + # Generated in Wolfram Language: + # + # Hyp[z_, p_ : 0, m_ : 0, s_ : 1, \[Lambda]_ : 1] := + # HyperbolicDistribution[\[Lambda], z*Sqrt[1 + p^2]/s, z*p/s, s, m] + # props[d_] := N[{Mean@d, Variance@d, Skewness@d, Kurtosis[d]-Kurtosis[NormalDistribution[]]}, 30] + # + # In[]:= props@Hyp[6, 1, -5, 2, 8] + + # Out[]= {1.06319517080168381227975950313, \ + # 5.63790040251113774016546423535, 0.625851630297306083508792293920, \ + # 0.63622204658230960216273474948} + distributions = [ + # No skewness, location, scale + ( + d=Hyp(3/10), + tmean=0, + tvar=23.7192375345642910012435489189, + tskew=0, + tkurt=2.68056465465234820166864415655, # EXCESS kurtosis!! + tmode=0 + ), ( + d=Hyp(3), + tmean=0, + tvar=0.510590348325343769505341346913, + tskew=0, + tkurt=0.88995291430008067418638784349, + tmode=0 + ), ( + d=Hyp(10), + tmean=0, + tvar=0.115341725074794522009169334735, + tskew=0, + tkurt=0.29539619644673402834977400441, + tmode=0 + ), + # Add skewness + ( + d=Hyp(1/10, -5), + tmean=-101.231534024878145278158003846, + tvar=10225.9463919953237264034118170, + tskew=-1.99460167233720176630548735540, + tkurt=5.97343971844269169149962193657, + tmode=-5 + ), ( + d=Hyp(3, -1), + tmean=-1.53177104497603130851602404074, + tvar=1.20662920739975591772951714814, + tskew=-1.14677758937397095739187195816, + tkurt=2.35384523161163514880079584352, + tmode=-1 + ), ( + d=Hyp(8, 1), + tmean=1.19272670576874667902532489770, + tvar=0.322857196451500261892719850967, + tskew=0.744755964647115580796052901480, + tkurt=1.05145730182810368503382261950, + tmode=1 + ), ( + d=Hyp(20, 5), + tmean=5.37946731533635375645988863585, + tvar=1.49459339171759999060447162252, + tskew=0.641583383199309715147192149775, + tkurt=0.68419269261165181164749648079, + tmode=5 + ), + # Add location & scale + ( + d=Hyp(1, -2, -1, 5), + tmean=-27.9948393559377234389267739977, + tvar=518.559320774469042942015734302, + tskew=-1.75030004538948425548612197692, + tkurt=4.77694919849117624423554369362, + tmode=-11 + ), ( + d=Hyp(1, -2, 1, 5), + tmean=-25.9948393559377234389267739977, + tvar=518.559320774469042942015734302, + tskew=-1.75030004538948425548612197692, + tkurt=4.77694919849117624423554369362, + tmode=-9 + ), ( + d=Hyp(6, 1, -5, 2), + tmean=-2.48204071280501662885409776329, + tvar=1.85647984002017020249703380834, + tskew=0.852297903911637048654316184661, + tkurt=1.35763888237097731740356582167, + tmode=-3 + ), + # Different λ + ( + d=Hyp(3/10, 0,0,1, -5), + tmean=0, + tvar=0.124766926919973221644770566113, + tskew=0, + tkurt=0.99265597462752447495230230680, + tmode=0 + ), ( + d=Hyp(3, 0,0,1, -1), + tmean=0, + tvar=0.288368126103121547283119124690, + tskew=0, + tkurt=1.00852166077397394811681109436, + tmode=0 + ), ( + d=Hyp(10, 0,0,1, 4), + tmean=0, + tvar=0.151980099806535291929364762119, + tskew=0, + tkurt=0.27275916387679590474273501914, + tmode=0 + ), ( + d=Hyp(3, -1, 0,1, -1), + tmean=-0.865104378309364641849357374071, + tvar=0.539962540733089251062850481473, + tskew=-1.21718592500810682678970904231, + tkurt=3.11624839625576195768292775000, + tmode=-0.560464077632561342794140059659 + ), ( + d=Hyp(8, 1, 0,1, 5), + tmean=1.84954386398016460385917258321, + tvar=0.584696274181089855715315261983, + tskew=0.668226904594935516683465663351, + tkurt=0.76850406887218305023668985068, + tmode=1.60212984628092067584739859051 + ), ( + d=Hyp(1, -2, -1, 5, -1/2), + tmean=-11, + tvar=125, + tskew=-2.68328157299974763569100840248, + tkurt=12.6, + tmode=-4.39143813989001345842141045503 + ), ( + d=Hyp(1, -2, -1, 5, 1/2), + tmean=-21, + tvar=350, + tskew=-2.02354940305121319832590492666, + tkurt=6.53571428571428571428571428571, + tmode=-7.49187863276682628498920382058 + ), ( + d=Hyp(6, 1, -5, 2, 8), + tmean=1.06319517080168381227975950313, + tvar=5.63790040251113774016546423535, + tskew=0.625851630297306083508792293920, + tkurt=0.63622204658230960216273474948, + tmode=0.327863171219295387089685071504 + ), + ] + NSAMPLES = 10^6 + for (d, mean_true, var_true, skew_true, kurt_true, mode_true) in distributions + println("\ttesting $d") + + @test collect(params(d)) ≈ [d.α, d.β, d.δ, d.μ, d.λ] + + # Specified `atol` because sometimes ground truth is zero, but we may get some floating-point inaccuracy + @test mean(d) ≈ mean_true atol=1e-6 + @test var(d) ≈ var_true atol=1e-6 + @test skewness(d) ≈ skew_true atol=1e-6 + @test kurtosis(d) ≈ kurt_true atol=1e-6 + @test mode(d) ≈ mode_true atol=1e-6 + + samples = rand(d, NSAMPLES) + @test abs(mean(d) - mean(samples)) < 0.02 + @test abs(std(d) - std(samples)) < 0.05 + + # Empirical CF should be close to theoretical CF + @test maximum(t->abs(cf(d, t) - cf_empirical(samples, t)), range(-50, 50, 100)) < 0.005 + + test_samples(d, NSAMPLES) + end +end diff --git a/test/univariate/continuous/generalizedinversegaussian.jl b/test/univariate/continuous/generalizedinversegaussian.jl new file mode 100644 index 000000000..350d45411 --- /dev/null +++ b/test/univariate/continuous/generalizedinversegaussian.jl @@ -0,0 +1,69 @@ +import SpecialFunctions: besselk + +@testset "Generalized inverse Gaussian" begin + GIG = GeneralizedInverseGaussian + # Empirical characteristic function + cf_empirical(samples::AbstractVector{<:Real}, t::Real) = mean(x->exp(1im * t * x), samples) + # Derivative d/dp log(besselk(p, x)) + dlog_besselk_dp(p::Real, x::Real, h::Real=1e-5) = (log(besselk(p + h, x)) - log(besselk(p - h, x))) / (2h) + + distributions = [ + # basic inverse Gaussian + ( + d=GIG(3, 8, -1/2), + tmean=3, + tvar=3.375, + tskew=1.83711730708738357364796305603, + tkurt=5.625 + ), ( + d=GIG(3, 8, 1.0), + tmean=4.80482214427053451026053858645, + tvar=7.53538381114490814886717269870, + tskew=1.46944720741022099816138306944, + tkurt=3.41236660328772156010103826417 + ), ( + d=GIG(1, 1, 1.0), + tmean=2.69948393559377234389267739977, + tvar=4.51072222384624734344698799307, + tskew=1.80431858973895795576206086081, + tkurt=4.99738802301388098591337070387 + ), ( + d=GIG(1, 1, -1.0), + tmean=0.699483935593772343892677399771, + tvar=0.510722223846247343446987993073, + tskew=3.52208925791368986307687496050, + tkurt=21.7757675389062712238907921144 + ) + ] + NSAMPLES = 10^6 + for (d, mean_true, var_true, skew_true, kurt_true) in distributions + println("\ttesting $d") + + @test collect(params(d)) ≈ [d.μ, d.λ, d.θ] + + @test mean(d) ≈ mean_true + @test var(d) ≈ var_true + @test skewness(d) ≈ skew_true + @test kurtosis(d) ≈ kurt_true + + samples = rand(d, NSAMPLES) + + @test abs(mean(d) - mean(samples)) < 0.01 + @test abs(std(d) - std(samples)) < 0.01 + + @test maximum(t->abs(cf(d, t) - cf_empirical(samples, t)), range(-50, 50, 100)) < 0.005 + + # a, b, p = params(d) + # t = sqrt(a * b) + # # E[log p(x; a,b,p)] + # expected_loglik_true = ( + # p/2 * log(a/b) - log(2besselk(p, t)) + # + (p-1) * (0.5 * log(b/a) + dlog_besselk_dp(p, t)) + # - (t * besselk(p+1, t)/besselk(p, t) - p) + # ) + # expected_loglik_sample = mean(x->logpdf(d, x), samples) + # @test abs(expected_loglik_true - expected_loglik_sample) < 0.01 + + test_samples(d, NSAMPLES) + end +end