Skip to content

Commit

Permalink
tests + cf + mgf
Browse files Browse the repository at this point in the history
  • Loading branch information
@jaksle
jaksle committed Dec 2, 2023
1 parent bb652b9 commit 768b727
Show file tree
Hide file tree
Showing 4 changed files with 154 additions and 13 deletions.
1 change: 1 addition & 0 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ version = "0.25.86"
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Expand Down
44 changes: 31 additions & 13 deletions src/univariate/continuous/stable.jl
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
# methods mainly from John P. Nolan, "Univariate Stable Distributions", Springer 2020

"""
Stable(α, β, σ, μ)
Expand Down Expand Up @@ -33,18 +34,18 @@ struct Stable{T<:Real} <: ContinuousUnivariateDistribution
β::T
σ::T
μ::T
Stable{T}(α, β, σ, μ) where {T} = new{T}(α, β, σ, μ)
function Stable{T}(α, β, σ, μ; check_args::Bool=true) where {T}
@check_args Stable (α, zero(T) < α <= 2one(T)) (β, -one(T) <= β <= one(T)) (σ, zero(T) < σ) ((α, β), α != 2 || β == 0 )
new{T}(α, β, σ, μ)
end
end

function Stable(α::T, β::T, σ::T, μ::T; check_args::Bool=true) where {T <: Real}
@check_args Stable (α, zero(T) < α <= 2one(T)) (β, -one(T) <= β <= one(T)) (σ, zero(T) < σ) ((α, β), α != 2 || β == 0 )
return Stable{T}(α, β, σ, μ)
end
Stable(α::T, β::T, σ::T, μ::T; check_args::Bool=true) where {T <: Real} = Stable{T}(α, β, σ, μ; check_args=check_args)

Stable(α::Real, β::Real, σ::Real, μ::Real; check_args::Bool=true) = Stable(promote(α, β, σ, μ)...; check_args=check_args)
Stable(α::Integer, β::Integer, σ::Integer, μ::Integer; check_args::Bool=true) = Stable(float(α), float(β), float(σ), float(μ); check_args=check_args)
Stable(α::Real, β::Real; check_args::Bool=true) = Stable(α, β, one(α), zero(α); check_args=check_args)
Stable(α::Real; check_args::Bool=true) = Stable(α,zero(α); check_args=check_args)
Stable(α::Real; check_args::Bool=true) = Stable(α, zero(α); check_args=check_args)

@distr_support Stable (d.α < 1 && d.β == 1 ? d.μ : -Inf) (d.α < 1 && d.β == -1 ? d.μ : Inf)

Expand All @@ -71,6 +72,23 @@ kurtosis(d::Stable{T}) where T = d.α == 2one(T) ? T(0.0) : T(NaN)

#### Evaluation

function cf(d::Stable{T}, t::Real) where T
α, β, σ, μ = params(d)
if α == one(T)
exp(im*t*μ - abs(σ*t) * (1 + im*β*2/π*sign(t)*log(abs(t))))
else
exp(im*t*μ - abs(σ*t)^α * (1 - im*β*sign(t)*tan(α*π/2)))
end
end

function mgf(d::Stable{T}, t::Real) where T
if d.α == 2one(T)
mgf(Normal(d.μ, √2d.σ), t)
else
T(Inf)
end
end

# integral representation from Nolan ch. 3
function pdf(d::Stable{T}, x::Real) where T
α, β, σ, μ = params(d)
Expand Down Expand Up @@ -126,7 +144,7 @@ function cdf(d::Stable{T}, x::Real) where T

z(v,c) = v*c > 36. ? 0.0 : exp(-c*v) # numerical truncation

function F(α,β,x) # works for x > 0
function F(α, β, x) # works for x > 0
V(θ) =(cos(α*θ₀))^(1/(α-1)) * (cos(θ)/sin(α*(θ₀+θ)))^(α/(α-1)) * cos(α*θ₀ + (α-1)*θ)/cos(θ)
θ₀= atan(β*tan(α*π/2))/α
x ≈ 0. && return (π/2 - θ₀)/π
Expand All @@ -136,20 +154,20 @@ function cdf(d::Stable{T}, x::Real) where T
return c + sign(1-α)/π * I
end

function F₁(β,x) # works for β > 0
function F₁(β, x) # works for β > 0
V₁(θ) = 2/π*(π/2+β*θ)/cos(θ) * exp((π/2+β*θ)*tan(θ)/β)
I, _err = quadgk(θ -> z(V₁(θ),exp(-π*x/2β)), -π/2, π/2 )
return 1/π* I
return 1/π * I
end

if α == one(T)
x = (x-μ)/σ - 2/π*β*log(σ) # normalize to S(1,β,1,0)
β < 0 && return F₁(-β,x)
return F₁(β,x)
β < 0 && return 1 - F₁(-β, -x)
return F₁(β, x)
else
x = (x-μ)/σ # normalize to S(α,β,1,0)
x < 0 && return 1 - F(α,-β,-x)
return F(α,β,x)
x < 0 && return 1 - F(α, -β, -x)
return F(α, β, x)

end
end
Expand Down
12 changes: 12 additions & 0 deletions test/convolution.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,18 @@ end
d4 = Gamma(1.2, 0.4)
@test_throws ArgumentError convolve(d1, d4)
end

@testset "Stable" begin
@test_throws ArgumentError convolve(Stable(1),Stable(1.5))
d1 = Stable(1.5, 1/2,2, 1)
d2 = Stable(1.5, -1/2, 1, 1)
d3 = @inferred(convolve(d1, d2))
@test d3 isa Stable
@test d3.α == 1.5
@test d3.β ≈ (d1.β*d1.σ^d3.α + d2.β*d2.σ^d3.α) / (d1.σ^d3.α + d2.σ^d3.α)
@test d3.σ ≈ (d1.σ^d3.α + d2.σ^d3.α)^(1/d3.α)
@test d3.μ == 1 + 1
end
end

@testset "continuous multivariate" begin
Expand Down
110 changes: 110 additions & 0 deletions test/univariate/continuous/stable.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,110 @@
using Test, Distributions

@testset "Stable" begin

@testset "input" begin
@test_throws DomainError Stable(-3.)
@test_throws DomainError Stable(3.)
@test_throws DomainError Stable(1., -2.)
@test_throws DomainError Stable(1., 2.)
@test_throws DomainError Stable(2, 0.1)
@test_throws DomainError Stable{Float64}(1.5, 0, -2, 0)

@test Stable(1, 1, 1, 1) isa Stable{Float64}
@test Stable(1.7) == Stable(1.7, 0, 1, 0)
@test Stable(0.2, -0.5) == Stable(0.2, -0.5, 1, 0)
end

@testset "conversions" begin
@test convert(Stable{Float64}, 1, 1, 1, 1) isa Stable{Float64}
@test convert(Stable{BigFloat}, Stable(1., 1, 1, 1)) isa Stable{BigFloat}
end

@testset "parameters" begin
@test shape(Stable(1.2, 1.)) == (1.2, 1.)
@test location(Stable(1., 0., 1., 3.)) == 3.
@test scale(Stable(1.5, 1., 2., 0.)) == 2.
@test params(Stable(0.3, 0.4, 0.5, 0.6)) == (0.3, 0.4, 0.5, 0.6)
@test partype(Stable{Int}(1, 1, 1, 1)) === Int
end

@testset "statistics and support" begin
@test mean(Stable(2., 0., 3., -5.)) == -5.
@test mean(Stable(1.1, 0., 3., 42.)) == 42.
@test isnan(mean(Stable(1, 1)))
@test var(Stable(2., 0., 3., 0.)) == 18.
@test isinf(var(Stable(1.9, 0., 3., 0.)))
@test iszero(skewness(Stable(2.0, 0., 2., 1.)))
@test isnan(skewness(Stable(1.5, 0., 2., 1.)))
@test iszero(kurtosis(Stable(2.0, 0., 2., 1.)))
@test isnan(kurtosis(Stable(1.5, 0., 2., 1.)))

@test support(Stable(1.1)) == RealInterval(-Inf, Inf)
@test support(Stable(1., 1.)) == RealInterval(-Inf, Inf)
@test support(Stable(0.9, 1.)) == RealInterval(0., Inf)
@test support(Stable(0.9, -1., 1., 2.)) == RealInterval(-Inf, 2.)
@test support(Stable(0.9, 0.9)) == RealInterval(-Inf, Inf)
@test support(Stable(0.9, -0.9, 1., 2.)) == RealInterval(-Inf, Inf)
end

# we cannot use test_affine_transformations at α == 1
@testset "affine transformations" begin
@test 2Stable(0.7, -1, 1, 1) + 1 == Stable(0.7, -1, 2, 3)
@test -2Stable(1.5, -1, 1, 1) + 1 == Stable(1.5, 1, 2, -1)
@test -3Stable(1, 1, 2, 1) + 2 ≈ Stable(1, -1, 6, -3 + 2 + 2/pi*6*log(3))
end

@testset "pdf, cdf and mgf in special cases" begin
xs = LinRange(-5, 5, 20)
@test pdf.(Stable(2., 0., 2., 1.), xs) ≈ pdf.(Normal(1., 2√2), xs)
@test cdf.(Stable(2., 0., 2., 1.), xs) ≈ cdf.(Normal(1., 2√2), xs)
@test pdf.(Stable(1., 0., 3., -1.), xs) == pdf.(Cauchy(-1., 3.), xs)
@test cdf.(Stable(1., 0., 3., -1.), xs) == cdf.(Cauchy(-1., 3.), xs)
@test pdf.(Stable(0.5, 1., 2., 1.), xs) == pdf.(Levy(1., 2.), xs)
@test cdf.(Stable(0.5, 1., 2., 1.), xs) == cdf.(Levy(1., 2.), xs)
@test pdf.(Stable(0.5, -1., 2., 3.), xs) == pdf.(Levy(-3., 2.), -xs)
@test cdf.(Stable(0.5, -1., 2., 3.), xs) ≈ 1 .- cdf.(Levy(-3., 2.), -xs)
@test mgf.(Stable(2., 0., 2., -2.), xs) == mgf.(Normal(-2., √2*2), xs)
end

# test values taken from Mathematica 13.2
@testset "pdf, cdf and cf in general case" begin
@test cf(Stable(0.5), -2) ≈ exp(-√2)
@test cf(Stable(1.5, -0.5, 0.5, 0.5), 4.2) ≈ exp(im*2.1 - (2.1)^1.5*(1 + im*0.5*tan(0.75*π)))
@test cf(Stable(1., -0.75, 1., 0.), -4.2) ≈ exp( -(4.2)*(1 + im*1.5/π*log(4.2)))
@test isinf(mgf(Stable(1.9), 42.))

# checking proper shape, Stable(α, β)
@test pdf(Stable(1.5), 0.) ≈ 0.287353 atol = 1e-6
@test pdf(Stable(1.1), 0.) ≈ 0.307141 atol = 1e-6
@test pdf(Stable(0.4), 0.) ≈ 1.057855 atol = 1e-6
@test pdf(Stable(0.2), 0.) ≈ 38.197186 atol = 1e-6
@test pdf(Stable(1.9, 1.), 3.) ≈ 0.0300991 atol = 1e-6
@test pdf(Stable(1.2, 1.), -1.) ≈ 0.0973176 atol = 1e-6
@test pdf(Stable(1., -0.9), -1.) ≈ 0.162803 atol = 1e-6
@test pdf(Stable(1., -0.5), -10.) ≈ 0.005098 atol = 1e-6
@test pdf(Stable(0.5, 0.5), -10.) ≈ 0.002181 atol = 1e-6
@test pdf(Stable(0.7, 1.), -0.01) ≈ 0.000000 atol = 1e-6

@test cdf(Stable(1.5), 0.) ≈ 0.500000 atol = 1e-6
@test cdf(Stable(0.2), 0.) ≈ 0.500000 atol = 1e-6
@test cdf(Stable(1.9, 1), 3.) ≈ 0.970814 atol = 1e-6
@test cdf(Stable(1.2, 1), -1.) ≈ 0.758715 atol = 1e-6
@test cdf(Stable(1, 0.5), -1.) ≈ 0.165443 atol = 1e-6
@test cdf(Stable(1., -0.9), -1.) ≈ 0.405070 atol = 1e-6
@test cdf(Stable(1., -0.5), -10.) ≈ 0.050327 atol = 1e-6
@test cdf(Stable(0.5, 0.5), -10.) ≈ 0.052669 atol = 1e-6
@test cdf(Stable(0.7, 1.), -0.01) ≈ 0.000000 atol = 1e-6

# checking affine transformations, Stable(α, β, σ, μ)
(σ, μ) = 3., -2.
@test pdf(Stable(1.1, 1., σ, μ), 2.5) ≈ 1/σ * pdf(Stable(1.1, 1), (2.5 - μ)/σ)
@test cdf(Stable(0.8, 0.1, σ, μ), 2.5) ≈ cdf(Stable(0.8, 0.1), (2.5 - μ)/σ)
@test pdf(Stable(1., -0.5, σ, μ), 2.5) ≈ 1/σ * pdf(Stable(1., -0.5), (2.5 - μ)/σ + 1/π*log(σ))
@test cdf(Stable(1., -0.5, σ, μ), 2.5) ≈ cdf(Stable(1., -0.5), (2.5 - μ)/σ + 1/π*log(σ))
@test cdf(Stable(1.4, 0.25, 2., -1.), 0.0) ≈ 0.697587 atol = 1e-6
@test cdf(Stable(1, 0.25, 2., -1.), 0.0) ≈ 0.581518 atol = 1e-6
@test pdf(Stable(1, 0.25, 2., -1.), 0.0) ≈ 0.128031 atol = 1e-6
end

end

0 comments on commit 768b727

Please sign in to comment.