From 67606ea4124a68e2b8ccbee8baaeef379ef1dccd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1s=20Lob=C3=A3o=20de=20Almeida?= Date: Wed, 4 Jan 2023 18:02:42 +0000 Subject: [PATCH] Correct alphastable rng for alpha in (1,2) --- src/AlphaStableDistributions.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/AlphaStableDistributions.jl b/src/AlphaStableDistributions.jl index abb42b2..20b834d 100644 --- a/src/AlphaStableDistributions.jl +++ b/src/AlphaStableDistributions.jl @@ -294,20 +294,28 @@ function Base.rand(rng::AbstractRNG, d::AlphaStable{T}) where {T<:AbstractFloat} α=d.α; β=d.β; scale=d.scale; loc=d.location (α < 0.1 || α > 2) && throw(DomainError(α, "α must be in the range 0.1 to 2")) abs(β) > 1 && throw(DomainError(β, "β must be in the range -1 to 1")) + ϕ = (rand(rng, T) - 0.5) * π + if α == one(T) && β == zero(T) return loc + scale * tan(ϕ) end - w = -log(rand(rng, T)) + + w = randexp(rng, T) + α == 2 && (return loc + 2*scale*sqrt(w)*sin(ϕ)) + β == zero(T) && (return loc + scale * ((cos((1-α)*ϕ) / w)^(one(T)/α - one(T)) * sin(α * ϕ) / cos(ϕ)^(one(T)/α))) + cosϕ = cos(ϕ) + if abs(α - one(T)) > 1e-8 ζ = β * tan(π * α / 2) aϕ = α * ϕ a1ϕ = (one(T) - α) * ϕ - return loc + scale * (( (sin(aϕ) + ζ * cos(aϕ))/cosϕ * ((cos(a1ϕ) + ζ*sin(a1ϕ))) / ((w*cosϕ)^((1-α)/α)) )) + return loc + scale * ((sin(aϕ) + ζ*cos(aϕ)) / (cosϕ^(1/α))) * ((cos(a1ϕ) + ζ*sin(a1ϕ)) / w) ^ ((one(T)-α)/α) end + bϕ = π/2 + β*ϕ x = 2/π * (bϕ * tan(ϕ) - β * log(π/2*w*cosϕ/bϕ)) α == one(T) || (x += β * tan(π*α/2))