Skip to content

Commit 61c767e

Browse files
ymtooymtoo
authored
Merge pull request #20 from ymtoo/alphastable-fit
Alphastable fit
2 parents 3b72837 + 6bbda71 commit 61c767e

File tree

3 files changed

+225
-43
lines changed

3 files changed

+225
-43
lines changed

‎Project.toml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ version = "1.0.0"
55

66
[deps]
77
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
8+
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
89
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
910
MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
1011
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

‎src/AlphaStableDistributions.jl‎

Lines changed: 194 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,9 @@ module AlphaStableDistributions
33
using LinearAlgebra, Statistics, Random
44
using StatsBase, Distributions, StaticArrays
55
using MAT, SpecialFunctions, ToeplitzMatrices
6+
using Interpolations
67

7-
8-
export AlphaStable, AlphaSubGaussian, fit
8+
export AlphaStable, SymmetricAlphaStable, AlphaSubGaussian, fit
99

1010
Base.@kwdef struct AlphaStable{T} <: Distributions.ContinuousUnivariateDistribution
1111
α::T = 1.5
@@ -33,45 +33,209 @@ Statistics.var(d::AlphaStable) = d.α == 2 ? 2d.scale^2 : Inf
3333
# mgf(d::AlphaStable, ::Any) = error("Not implemented")
3434
# cf(d::AlphaStable, ::Any) = error("Not implemented")
3535

36-
# lookup table from McCulloch (1986)
36+
# lookup tables from McCulloch (1986)
3737
const _ena = [
38-
2.4388
39-
2.5120
40-
2.6080
41-
2.7369
42-
2.9115
43-
3.1480
44-
3.4635
45-
3.8824
46-
4.4468
47-
5.2172
48-
6.3140
49-
7.9098
50-
10.4480
51-
14.8378
52-
23.4831
53-
44.2813
38+
2.4388
39+
2.5120
40+
2.6080
41+
2.7369
42+
2.9115
43+
3.1480
44+
3.4635
45+
3.8824
46+
4.4468
47+
5.2172
48+
6.3140
49+
7.9098
50+
10.4480
51+
14.8378
52+
23.4831
53+
44.2813
54+
]
55+
const _να = [
56+
2.439
57+
2.500
58+
2.600
59+
2.700
60+
2.800
61+
3.000
62+
3.200
63+
3.500
64+
4.000
65+
5.000
66+
6.000
67+
8.000
68+
10.000
69+
15.000
70+
25.000
71+
]
72+
const _νβ = [
73+
0.0
74+
0.1
75+
0.2
76+
0.3
77+
0.5
78+
0.7
79+
1.0
80+
]
81+
const ψ₁ = [
82+
2.000 2.000 2.000 2.000 2.000 2.000 2.000
83+
1.916 1.924 1.924 1.924 1.924 1.924 1.924
84+
1.808 1.813 1.829 1.829 1.829 1.829 1.829
85+
1.729 1.730 1.737 1.745 1.745 1.745 1.745
86+
1.664 1.663 1.663 1.668 1.676 1.676 1.676
87+
1.563 1.560 1.553 1.548 1.547 1.547 1.547
88+
1.484 1.480 1.471 1.460 1.448 1.438 1.438
89+
1.391 1.386 1.378 1.364 1.337 1.318 1.318
90+
1.279 1.273 1.266 1.250 1.210 1.184 1.150
91+
1.128 1.121 1.114 1.101 1.067 1.027 0.973
92+
1.029 1.021 1.014 1.004 0.974 0.935 0.874
93+
0.896 0.892 0.887 0.883 0.855 0.823 0.769
94+
0.818 0.812 0.806 0.801 0.780 0.756 0.691
95+
0.698 0.695 0.692 0.689 0.676 0.656 0.595
96+
0.593 0.590 0.588 0.586 0.579 0.563 0.513
97+
]
98+
const ψ₂ = [
99+
0.0 2.160 1.000 1.000 1.000 1.000 1.000
100+
0.0 1.592 3.390 1.000 1.000 1.000 1.000
101+
0.0 0.759 1.800 1.000 1.000 1.000 1.000
102+
0.0 0.482 1.048 1.694 2.229 1.000 1.000
103+
0.0 0.360 0.760 1.232 2.229 1.000 1.000
104+
0.0 0.253 0.518 0.823 1.575 1.000 1.000
105+
0.0 0.203 0.410 0.632 1.244 1.906 1.000
106+
0.0 0.165 0.332 0.499 0.943 1.560 1.000
107+
0.0 0.136 0.271 0.404 0.689 1.230 2.195
108+
0.0 0.109 0.216 0.323 0.539 0.827 1.917
109+
0.0 0.096 0.190 0.284 0.472 0.693 1.759
110+
0.0 0.082 0.163 0.243 0.412 0.601 1.596
111+
0.0 0.074 0.147 0.220 0.377 0.546 1.482
112+
0.0 0.064 0.128 0.191 0.330 0.478 1.362
113+
0.0 0.056 0.112 0.167 0.285 0.428 1.274
114+
]
115+
const _α =[
116+
0.5
117+
0.6
118+
0.7
119+
0.8
120+
0.9
121+
1.0
122+
1.1
123+
1.2
124+
1.3
125+
1.4
126+
1.5
127+
1.6
128+
1.7
129+
1.8
130+
1.9
131+
2.0
132+
]
133+
const _β = [
134+
0.0
135+
0.25
136+
0.50
137+
0.75
138+
1.00
139+
]
140+
const ϕ₃ = [
141+
2.588 3.073 4.534 6.636 9.144
142+
2.337 2.635 3.542 4.808 6.247
143+
2.189 2.392 3.004 3.844 4.775
144+
2.098 2.244 2.676 3.265 3.912
145+
2.040 2.149 2.461 2.886 3.356
146+
2.000 2.085 2.311 2.624 2.973
147+
1.980 2.040 2.205 2.435 2.696
148+
1.965 2.007 2.125 2.294 2.491
149+
1.955 1.984 2.067 2.188 2.333
150+
1.946 1.967 2.022 2.106 2.211
151+
1.939 1.952 1.988 2.045 2.116
152+
1.933 1.940 1.962 1.997 2.043
153+
1.927 1.930 1.943 1.961 1.987
154+
1.921 1.922 1.927 1.936 1.947
155+
1.914 1.915 1.916 1.918 1.921
156+
1.908 1.908 1.908 1.908 1.908
157+
]
158+
const ϕ₅ = [
159+
0.0 0.0 0.0 0.0 0.0
160+
0.0 -0.017 -0.032 -0.049 -0.064
161+
0.0 -0.030 -0.061 -0.092 -0.123
162+
0.0 -0.043 -0.088 -0.132 -0.179
163+
0.0 -0.056 -0.111 -0.170 -0.232
164+
0.0 -0.066 -0.134 -0.206 -0.283
165+
0.0 -0.075 -0.154 -0.241 -0.335
166+
0.0 -0.084 -0.173 -0.276 -0.390
167+
0.0 -0.090 -0.192 -0.310 -0.447
168+
0.0 -0.095 -0.208 -0.346 -0.508
169+
0.0 -0.098 -0.223 -0.383 -0.576
170+
0.0 -0.099 -0.237 -0.424 -0.652
171+
0.0 -0.096 -0.250 -0.469 -0.742
172+
0.0 -0.089 -0.262 -0.520 -0.853
173+
0.0 -0.078 -0.272 -0.581 -0.997
174+
0.0 -0.061 -0.279 -0.659 -1.198
54175
]
55176

56177
"""
57178
fit(d::Type{<:AlphaStable}, x; alg=QuickSort)
58179
59-
Fit a symmetric α stable distribution to data.
180+
Fit an α stable distribution to data.
60181
61182
returns `AlphaStable`
62183
63-
α, c and δ are the characteristic exponent, scale parameter
184+
α∈[0.6,2.0], β∈[-1,1] , c∈[0,∞] and δ∈[-∞,∞] are the characteristic exponent,
185+
skewness parameter, scale parameter (dispersion^1/α) and location parameter respectively.
186+
187+
α, β, c and δ are computed based on McCulloch (1986) fractile.
188+
"""
189+
function Distributions.fit(::Type{<:AlphaStable}, x)
190+
p = quantile.(Ref(sort(x)), (0.05, 0.25, 0.28, 0.5, 0.72, 0.75, 0.95), sorted=true)
191+
να = (p[7]-p[1]) / (p[6]-p[2])
192+
νβ = (p[7]+p[1]-2p[4]) / (p[7]-p[1])
193+
(να < _να[1]) && (να = _να[1])
194+
(να > _να[end]) && (να = _να[end])
195+
196+
itp₁ = interpolate((_να, _νβ), ψ₁, Gridded(Linear()))
197+
α = itp₁(να, abs(νβ))
198+
itp₂ = interpolate((_να, _νβ), ψ₂, Gridded(Linear()))
199+
β = sign(νβ) * itp₂(να, abs(νβ))
200+
201+
(β > 1.0) && (β = 1.0)
202+
(β < -1.0) && (β = -1.0)
203+
itp₃ = interpolate((_α, _β), ϕ₃, Gridded(Linear()))
204+
c = (p[6]-p[2]) / itp₃(α, abs(β))
205+
itp₄ = interpolate((_α, _β), ϕ₅, Gridded(Linear()))
206+
ζ = p[4] + c * sign(β) * itp₄(α, abs(β))
207+
if abs(α - 1.0) < 0.1
208+
δ = ζ
209+
else
210+
δ = ζ - β * c * tan(π*α/2)
211+
end
212+
return AlphaStable(α=α, β=β, scale=c, location=oftype(α, δ))
213+
end
214+
215+
Base.@kwdef struct SymmetricAlphaStable{T} <: Distributions.ContinuousUnivariateDistribution
216+
α::T = 1.5
217+
scale::T = one(α)
218+
location::T = zero(α)
219+
end
220+
221+
"""
222+
fit(d::Type{<:SymmetricAlphaStable}, x; alg=QuickSort)
223+
224+
Fit a symmetric α stable distribution to data.
225+
226+
returns `SymmetricAlphaStable`
227+
228+
α∈[1,2], c∈[0,∞] and δ∈[-∞,∞] are the characteristic exponent, scale parameter
64229
(dispersion^1/α) and location parameter respectively.
65230
66231
α is computed based on McCulloch (1986) fractile.
67-
c is computed based on Fama & Roll (1971) fractile.
68-
δ is the 50% trimmed mean of the sample.
232+
scale is computed based on Fama & Roll (1971) fractile.
233+
location is the 50% trimmed mean of the sample.
69234
"""
70-
function Distributions.fit(d::Type{<:AlphaStable}, x; alg=QuickSort)
71-
sx = sort(x, alg=alg)
72-
δ = mean(@view(sx[end÷4:(3*end)÷4]))
73-
p = quantile.(Ref(sx), (0.05, 0.25, 0.28, 0.72, 0.75, 0.95), sorted=true)
74-
c = (p[4]-p[3])/1.654
235+
function Distributions.fit(::Type{<:SymmetricAlphaStable}, x)
236+
δ = mean(StatsBase.trim(x,prop=0.25))
237+
p = quantile.(Ref(sort(x)), (0.05, 0.25, 0.28, 0.72, 0.75, 0.95), sorted=true)
238+
c = (p[4]-p[3])/1.654
75239
an = (p[6]-p[1])/(p[5]-p[2])
76240
if an < 2.4388
77241
α = 2.
@@ -86,7 +250,7 @@ function Distributions.fit(d::Type{<:AlphaStable}, x; alg=QuickSort)
86250
if α < 0.5
87251
α = 0.5
88252
end
89-
return AlphaStable(α=α, β=zero(α), scale=c, location=oftype(α, δ))
253+
return SymmetricAlphaStable(α=α, scale=c, location=oftype(α, δ))
90254
end
91255

92256
"""
@@ -295,4 +459,4 @@ function Distributions.fit(d::Type{<:AlphaSubGaussian}, x::AbstractVector{T}, m:
295459
AlphaSubGaussian(α=α, R=cov, n=length(x))
296460
end
297461

298-
end # module
462+
end # module

‎test/runtests.jl‎

Lines changed: 30 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,22 +1,39 @@
11
using AlphaStableDistributions
2-
using Test, Random
2+
using Test, Random, Distributions
33

44
@testset "AlphaStableDistributions.jl" begin
55

6-
for α in 0.5:0.1:2
7-
d1 = AlphaStable(α=α)
8-
s = rand(d1, 100000)
6+
stabletypes = [AlphaStable,SymmetricAlphaStable]
7+
αs = [0.6:0.1:2,1:0.1:2]
8+
for (i, stabletype) in enumerate(stabletypes)
9+
for α in αs[i]
10+
d1 = AlphaStable(α=α)
11+
s = rand(d1, 100000)
912

10-
d2 = fit(AlphaStable, s)
13+
d2 = fit(stabletype, s)
1114

12-
@test d1.α ≈ d2.α rtol=0.1
13-
@test d1.β ≈ d2.β rtol=0.1
14-
@test d1.scale ≈ d2.scale rtol=0.1
15-
@test d1.location ≈ d2.location atol=0.03
15+
@test d1.α ≈ d2.α rtol=0.1
16+
stabletype != SymmetricAlphaStable && @test d1.β ≈ d2.β atol=0.2
17+
@test d1.scale ≈ d2.scale rtol=0.1
18+
@test d1.location ≈ d2.location atol=0.1
19+
end
20+
21+
xnormal = rand(Normal(3.0, 4.0), 96000)
22+
d = fit(stabletype, xnormal)
23+
@test d.α ≈ 2 rtol=0.2
24+
stabletype != SymmetricAlphaStable && @test d.β ≈ 0 atol=0.2
25+
@test d.scale ≈ 4/√2 rtol=0.2
26+
@test d.location ≈ 3 rtol=0.1
27+
28+
xcauchy = rand(Cauchy(3.0, 4.0), 96000)
29+
d = fit(stabletype, xcauchy)
30+
@test d.α ≈ 1 rtol=0.2
31+
stabletype != SymmetricAlphaStable && @test d.β ≈ 0 atol=0.2
32+
@test d.scale ≈ 4 rtol=0.2
33+
@test d.location ≈ 3 rtol=0.1
1634
end
1735

1836
for α in 1.1:0.1:1.9
19-
2037
d = AlphaSubGaussian(n=96000, α=α)
2138
x = rand(d)
2239
x2 = copy(x)
@@ -25,9 +42,9 @@ using Test, Random
2542

2643
d3 = fit(AlphaStable, x)
2744
@test d3.α ≈ α rtol=0.2
28-
@test d3.β == 0
45+
@test d3.β ≈ 0 atol=0.2
2946
@test d3.scale ≈ 1 rtol=0.2
30-
@test d3.location ≈ 0 atol=0.03
47+
@test d3.location ≈ 0 atol=0.1
3148
end
3249

3350
d4 = AlphaSubGaussian(n=96000)
@@ -91,4 +108,4 @@ end
91108
# d1 = AlphaStable(α=1.5)
92109
# s = rand(d1, 100000)
93110
# using ThreadsX
94-
# @btime fit($AlphaStable, $s, $ThreadsX.MergeSort)
111+
# @btime fit($AlphaStable, $s, $ThreadsX.MergeSort)

0 commit comments

Comments
 (0)