Merge pull request #35 from jw3126/cf

implement cf for (Symmetric)AlphaStable

Author: ymtoo

src/AlphaStableDistributions.jl

@@ -15,7 +15,11 @@ Base.@kwdef struct AlphaStable{T} <: Distributions.ContinuousUnivariateDistribut
 end
 
 AlphaStable(α::Integer, β::Integer, scale::Integer, location::Integer) = AlphaStable(float(α), float(β), float(scale), float(location))
-
+function AlphaStable(α,β,scale,location)
+    αT,βT,scaleT,locationT =  promote(α,β,scale,location)
+    AlphaStable(αT,βT,scaleT,locationT)
+end
+Distributions.params(d::AlphaStable) = (d.α, d.β, d.scale, d.location)
 
 # sampler(d::AlphaStable) = error("Not implemented")
 # pdf(d::AlphaStable, x::Real) = error("Not implemented")
@@ -33,7 +37,17 @@ Statistics.var(d::AlphaStable) = d.α == 2 ? 2d.scale^2 : Inf
 # kurtosis(d::Distribution, ::Bool) = error("Not implemented")
 # entropy(d::AlphaStable, ::Real) = error("Not implemented")
 # mgf(d::AlphaStable, ::Any) = error("Not implemented")
-# cf(d::AlphaStable, ::Any) = error("Not implemented")
+
+function Distributions.cf(d::AlphaStable{S}, t::Real) where {S}
+    T = float(promote_type(S, typeof(t)))
+    α,β,c,δ = params(d)
+    Φ = if α == one(α)
+        T(-2/π * log(abs(t)))
+    else
+        T(tan(π*α/2))
+    end
+    exp(im*t*δ - abs(c*t)^α * (1 - im*β*sign(t)*Φ))
+end
 
 # lookup tables from McCulloch (1986)
 const _ena = [
@@ -220,6 +234,12 @@ Base.@kwdef struct SymmetricAlphaStable{T} <: Distributions.ContinuousUnivariate
     scale::T = one(α)
     location::T = zero(α)
 end
+Distributions.params(d::SymmetricAlphaStable) = (d.α, d.scale, d.location)
+Distributions.cf(d::SymmetricAlphaStable, t::Real) = cf(AlphaStable(d), t)
+
+function AlphaStable(d::SymmetricAlphaStable)
+    AlphaStable(α=d.α,scale=d.scale,location=d.location)
+end
 
 """
     fit(d::Type{<:SymmetricAlphaStable}, x; alg=QuickSort)
@@ -465,4 +485,4 @@ function Distributions.fit(d::Type{<:AlphaSubGaussian}, x::AbstractVector{T}, m:
     AlphaSubGaussian(α=α, R=cov, n=length(x))
 end
 
-end # module
+end # module

test/runtests.jl

@@ -1,16 +1,80 @@
 using AlphaStableDistributions
 using Test, Random, Distributions
 
-@testset "AlphaStableDistributions.jl" begin
+@testset "cf" begin
+    rng = MersenneTwister(1)
+    for _ in 1:100
+        d = AlphaStable(
+            α=rand(rng,Uniform(0,2)), 
+            β=rand(rng,Uniform(-1,1)), 
+            scale=rand(rng,Uniform(0,10)), 
+            location=rand(rng,Uniform(0,10)),
+        )
+        @test cf(d,0) ≈ 1
+        x = rand(rng,Uniform(-10,10))
+        @test abs(cf(d, x)) <= 1
+
+        d32 = AlphaStable(Float32.(Distributions.params(d))...)
+        @test cf(d32, Float32(x)) isa Complex{Float32}
+        @test cf(d32, Float32(x)) ≈ cf(d,x) atol=100*eps(Float32)
+    end
+    for _ in 1:100
+        # test stability under convolution
+        d = AlphaStable(
+            α=rand(rng,Uniform(0.1,2)), 
+            scale=1.0,
+            location=0.0,
+        )
+        x = rand(rng,Uniform(-1,1))
+        n = rand(rng,1:10^6)
+        s = n^inv(d.α)
+        @test cf(d, x) ≈ cf(d, x/s)^n
+    end
+    xs = range(-1,1,length=100)
+    d1 = SymmetricAlphaStable(α=2.0, scale=1/sqrt(2), location=0.0)
+    d2 = AlphaStable(d1)
+    d_ref = Normal(0.0,1.0)
+    @test cf.(Ref(d1),xs) ≈ cf.(Ref(d_ref),xs)
+    @test cf.(Ref(d2),xs) ≈ cf.(Ref(d_ref),xs)
+    
+    xs = range(-10,10,length=100)
+    d1 = SymmetricAlphaStable(α=1.0, scale=1.0, location=0.0)
+    d2 = AlphaStable(d1)
+    d_ref = Cauchy(0.0,1.0)
+    @test cf.(Ref(d1),xs) ≈ cf.(Ref(d_ref),xs)
+    @test cf.(Ref(d2),xs) ≈ cf.(Ref(d_ref),xs)
+    
+    d1 = SymmetricAlphaStable(α=1.0, scale=17.9, location=42.0)
+    d2 = AlphaStable(d1)
+    d_ref = Cauchy(42.0,17.9)
+    @test cf.(Ref(d1),xs) ≈ cf.(Ref(d_ref),xs)
+    @test cf.(Ref(d2),xs) ≈ cf.(Ref(d_ref),xs)
+    
+    d1 = AlphaStable(α=1/2, β=1.0, scale=12.0, location=-7.2)
+    d_ref = Levy(-7.2, 12.0)
+    @test cf.(Ref(d1),xs) ≈ cf.(Ref(d_ref),xs)
 
+    @test @inferred(cf(AlphaStable(α = 1.0) , 1.0)) isa Complex{Float64}
+    @test @inferred(cf(AlphaStable(α = 1  ) , 1  )) isa Complex{Float64}
+    @test @inferred(cf(AlphaStable(α = 1.0) , 1f0)) isa Complex{Float64}
+    @test @inferred(cf(AlphaStable(α = 1f0) , 1  )) isa Complex{Float32}
+    @test @inferred(cf(AlphaStable(α = 1f0) , 1f0)) isa Complex{Float32}
+end
+
+    
+@testset "AlphaStableDistributions.jl" begin
+    @test AlphaStable(α=1, scale=1.5) === AlphaStable(α=1.0, scale=1.5)
+    @test Distributions.params(AlphaStable()) === (1.5, 0.0, 1.0, 0.0)
+    @test Distributions.params(SymmetricAlphaStable()) === (1.5, 1.0, 0.0)
+    rng = MersenneTwister(0)
     sampletypes = [Float32,Float64]
     stabletypes = [AlphaStable,SymmetricAlphaStable]
     αs = [0.6:0.1:2,1:0.1:2]
     for sampletype ∈ sampletypes
         for (i, stabletype) in enumerate(stabletypes)
             for α in αs[i]
                 d1 = AlphaStable(α=sampletype(α))
-                s = rand(d1, 200000)
+                s = rand(rng,d1, 200000)
                 @test eltype(s) == sampletype
 
                 d2 = fit(stabletype, s)
@@ -22,14 +86,14 @@ using Test, Random, Distributions
                 @test d1.location ≈ d2.location atol=0.1
             end
 
-            xnormal = rand(Normal(3.0, 4.0), 96000)
+            xnormal = rand(rng,Normal(3.0, 4.0), 96000)
             d = fit(stabletype, xnormal)
             @test d.α ≈ 2 rtol=0.2
             stabletype != SymmetricAlphaStable && @test d.β ≈ 0 atol=0.2
             @test d.scale ≈ 4/√2 rtol=0.2
             @test d.location ≈ 3 rtol=0.1
 
-            xcauchy = rand(Cauchy(3.0, 4.0), 96000)
+            xcauchy = rand(rng,Cauchy(3.0, 4.0), 96000)
             d = fit(stabletype, xcauchy)
             @test d.α ≈ 1 rtol=0.2
             stabletype != SymmetricAlphaStable && @test d.β ≈ 0 atol=0.2
@@ -40,9 +104,9 @@ using Test, Random, Distributions
 
     for α in 1.1:0.1:1.9
         d = AlphaSubGaussian(α=α, n=96000)
-        x = rand(d)
+        x = rand(rng,d)
         x2 = copy(x)
-        rand!(d, x2)
+        rand!(rng,d, x2)
         @test x != x2
 
         d3 = fit(AlphaStable, x)
@@ -54,13 +118,12 @@ using Test, Random, Distributions
 
     d4 = AlphaSubGaussian(α=1.5, n=96000)
     m = size(d4.R, 1) - 1
-    x = rand(d4)
+    x = rand(rng,d4)
     d5 = fit(AlphaSubGaussian, x, m, p=1.0)
     @test d4.α ≈ d5.α rtol=0.1
     @test d4.R ≈ d5.R rtol=0.1
 
 end
-
 # 362.499 ms (4620903 allocations: 227.64 MiB)
 # 346.520 ms (4621052 allocations: 209.62 MiB) # StaticArrays in outer fun
 # 345.925 ms (4225524 allocations: 167.66 MiB) # tempind to tuple
@@ -114,4 +177,4 @@ end
 # d1 = AlphaStable(α=1.5)
 # s = rand(d1, 100000)
 # using ThreadsX
-# @btime fit($AlphaStable, $s, $ThreadsX.MergeSort)
+# @btime fit($AlphaStable, $s, $ThreadsX.MergeSort)