diff --git a/Project.toml b/Project.toml index 0e897b4..b07d245 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "AlphaStableDistributions" uuid = "f20549b4-2d50-407f-863c-cdd202ba59a3" authors = ["Fredrik Bagge Carlson", "Too Yuen Min"] -version = "0.1.2" +version = "1.0.0" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/README.md b/README.md index 5e56789..a606bec 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,7 @@ AlphaStable{Float64}(α=1.5, β=0.0, scale=1.0, location=0.0) julia> s = rand(d1, 100000); -julia> d2 = fit(AlphaStable, s) +julia> d2 = fit(AlphaStable, s, alg=QuickSort) # See ThreadsX.QuickSort for a threaded algorithm AlphaStable{Float64}(α=1.4748701622930906, β=0.0, scale=1.006340087707924, location=-0.0036724481641865715) diff --git a/src/AlphaStableDistributions.jl b/src/AlphaStableDistributions.jl index 832e14c..26c75ae 100644 --- a/src/AlphaStableDistributions.jl +++ b/src/AlphaStableDistributions.jl @@ -54,10 +54,11 @@ const _ena = [ ] """ + fit(d::Type{<:AlphaStable}, x; alg=QuickSort) + Fit a symmetric α stable distribution to data. -:param x: data -:returns: (α, c, δ) +returns `AlphaStable` α, c and δ are the characteristic exponent, scale parameter (dispersion^1/α) and location parameter respectively. @@ -66,11 +67,11 @@ Fit a symmetric α stable distribution to data. c is computed based on Fama & Roll (1971) fractile. δ is the 50% trimmed mean of the sample. """ -function Distributions.fit(d::Type{<:AlphaStable}, x, alg=QuickSort) +function Distributions.fit(d::Type{<:AlphaStable}, x; alg=QuickSort) sx = sort(x, alg=alg) - δ = mean(@view(sx[end÷4:(3*end)÷4])) - p = quantile.(Ref(sx), (0.05, 0.25, 0.28, 0.72, 0.75, 0.95), sorted=true) - c = (p[4]-p[3])/1.654 + δ = mean(@view(sx[end÷4:(3*end)÷4])) + p = quantile.(Ref(sx), (0.05, 0.25, 0.28, 0.72, 0.75, 0.95), sorted=true) + c = (p[4]-p[3])/1.654 an = (p[6]-p[1])/(p[5]-p[2]) if an < 2.4388 α = 2. diff --git a/test/runtests.jl b/test/runtests.jl index d655401..0ecf0f2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -86,3 +86,9 @@ end # using Cthulhu # d = AlphaSubGaussian(n=96) # @descend_code_warntype rand(Random.GLOBAL_RNG, d) +# +# using AlphaStableDistributions +# d1 = AlphaStable(α=1.5) +# s = rand(d1, 100000) +# using ThreadsX +# @btime fit($AlphaStable, $s, $ThreadsX.MergeSort)