Merge pull request #47 from org-arl/fix-rand

Fix bugs

Author: ymtoo

Project.toml

@@ -1,7 +1,7 @@
 name = "AlphaStableDistributions"
 uuid = "f20549b4-2d50-407f-863c-cdd202ba59a3"
 authors = ["Fredrik Bagge Carlson", "Too Yuen Min"]
-version = "1.1.5"
+version = "1.1.6"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/AlphaStableDistributions.jl

@@ -172,22 +172,22 @@ const ϕ₃ = [
     1.908  1.908  1.908  1.908  1.908
 ]
 const ϕ₅ = [
-    0.0 0.0    0.0    0.0    0.0
-    0.0 -0.017 -0.032 -0.049 -0.064
-    0.0 -0.030 -0.061 -0.092 -0.123
-    0.0 -0.043 -0.088 -0.132 -0.179
-    0.0 -0.056 -0.111 -0.170 -0.232
-    0.0 -0.066 -0.134 -0.206 -0.283
-    0.0 -0.075 -0.154 -0.241 -0.335
-    0.0 -0.084 -0.173 -0.276 -0.390
-    0.0 -0.090 -0.192 -0.310 -0.447
-    0.0 -0.095 -0.208 -0.346 -0.508
-    0.0 -0.098 -0.223 -0.383 -0.576
-    0.0 -0.099 -0.237 -0.424 -0.652
-    0.0 -0.096 -0.250 -0.469 -0.742
-    0.0 -0.089 -0.262 -0.520 -0.853
-    0.0 -0.078 -0.272 -0.581 -0.997
-    0.0 -0.061 -0.279 -0.659 -1.198
+    0.0  -0.061  -0.279  -0.659  -1.198
+    0.0  -0.078  -0.272  -0.581  -0.997
+    0.0  -0.089  -0.262  -0.52   -0.853
+    0.0  -0.096  -0.25   -0.469  -0.742
+    0.0  -0.099  -0.237  -0.424  -0.652
+    0.0  -0.098  -0.223  -0.383  -0.576
+    0.0  -0.095  -0.208  -0.346  -0.508
+    0.0  -0.09   -0.192  -0.31   -0.447
+    0.0  -0.084  -0.173  -0.276  -0.39
+    0.0  -0.075  -0.154  -0.241  -0.335
+    0.0  -0.066  -0.134  -0.206  -0.283
+    0.0  -0.056  -0.111  -0.17   -0.232
+    0.0  -0.043  -0.088  -0.132  -0.179
+    0.0  -0.03   -0.061  -0.092  -0.123
+    0.0  -0.017  -0.032  -0.049  -0.064
+    0.0   0.0     0.0     0.0     0.0
 ]
 
 """
@@ -221,7 +221,7 @@ function Distributions.fit(::Type{<:AlphaStable}, x::AbstractArray{T}, alg=Quick
     c = (p[6]-p[2]) / itp₃(α, abs(β))
     itp₄ = interpolate((_α, _β), ϕ₅, Gridded(Linear()))
     ζ = p[4] + c * sign(β) * itp₄(α, abs(β))
-    if abs(α - 1.0) < 0.1
+    if abs(α - 1.0) < 0.05
         δ = ζ
     else
         δ = ζ - β * c * tan(π*α/2)
@@ -237,6 +237,7 @@ end
 Distributions.params(d::SymmetricAlphaStable) = (d.α, d.scale, d.location)
 Distributions.cf(d::SymmetricAlphaStable, t::Real) = cf(AlphaStable(d), t)
 Random.rand(rng::AbstractRNG, d::SymmetricAlphaStable) = rand(rng, AlphaStable(d))
+Base.eltype(::Type{<:SymmetricAlphaStable{T}}) where {T<:AbstractFloat} = T
 
 function AlphaStable(d::SymmetricAlphaStable)
     AlphaStable(α=d.α,scale=d.scale,location=d.location)
@@ -290,29 +291,29 @@ This implementation is based on the method in J.M. Chambers, C.L. Mallows
 and B.W. Stuck, "A Method for Simulating Stable Random Variables," JASA 71 (1976): 340-4.
 McCulloch's MATLAB implementation (1996) served as a reference in developing this code.
 """
-function Base.rand(rng::AbstractRNG, d::AlphaStable{T}) where {T<:AbstractFloat}
-    α=d.α; β=d.β; scale=d.scale; loc=d.location
+function Base.rand(rng::AbstractRNG, d::AlphaStable{T}) where {T<:AbstractFloat} 
+    α=d.α; β=d.β; sc=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) * π
+    # added eps(T) to prevent DomainError: x ^ y where x < 0
+    ϕ = (rand(rng, T) - T(0.5)) * π * (one(T) - eps(T))
     if α == one(T) && β == zero(T)
-        return loc + scale * tan(ϕ)
+        return loc + sc * tan(ϕ)
     end
     w = -log(rand(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)/α)))
+    α == 2 && (return loc + 2*sc*sqrt(w)*sin(ϕ))
+    β == zero(T) && (return loc + sc * ((cos((one(T)-α)*ϕ) / 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 + sc * ((((sin(aϕ) + ζ * cos(aϕ))/cosϕ) * ((cos(a1ϕ) + ζ*sin(a1ϕ)) / (w*cosϕ))^((one(T)-α)/α)))
     end
     bϕ = π/2 + β*ϕ
     x = 2/π * (bϕ * tan(ϕ) - β * log(π/2*w*cosϕ/bϕ))
     α == one(T) || (x += β * tan(π*α/2))
-
-    return loc + scale * x
+    return loc + sc * x
 end
 
 Base.eltype(::Type{<:AlphaStable{T}}) where {T<:AbstractFloat} = T
@@ -466,10 +467,10 @@ with memory", Signal Processing, Volume 131, Pages 271-279, 2017.
 """
 function Distributions.fit(d::Type{<:AlphaSubGaussian}, x::AbstractVector{T}, m::Integer; p=one(T)) where T
     d1   = fit(AlphaStable, x)
-    α    = d1.α; scale=d1.scale
+    α    = d1.α; sc=d1.scale
     cov  = zeros(T, m+1, m+1)
     xlen = length(x)
-    c    = ((sum(x->abs(x)^p, x)/xlen)^(1/p))/scale
+    c    = ((sum(x->abs(x)^p, x)/xlen)^(1/p))/sc
     for i in 1:m
         tempxlen = xlen-mod(xlen, i)
         xtemp = reshape(x[1:end-mod(xlen, i)], i, tempxlen÷i)
@@ -478,11 +479,11 @@ function Distributions.fit(d::Type{<:AlphaSubGaussian}, x::AbstractVector{T}, m:
             tempxlen = size(xtemp, 1)*size(xtemp, 2)
         end
         xtemp = reshape(xtemp', 2, tempxlen÷2)
-        @views r = (2/(c^p))*(scale^(2-p))*(xtemp[1, :]'*((sign.(xtemp[2, :]).*(abs.(xtemp[2, :]).^(p-1)))))/(tempxlen/2)
+        @views r = (2/(c^p))*(sc^(2-p))*(xtemp[1, :]'*((sign.(xtemp[2, :]).*(abs.(xtemp[2, :]).^(p-1)))))/(tempxlen/2)
         cov[diagind(cov, i)] .+= r
     end
-    cov = (cov+cov')+2*(scale^2)*I(m+1)
-    cov ./= 2*scale^2
+    cov = (cov+cov')+2*(sc^2)*I(m+1)
+    cov ./= 2*sc^2
     AlphaSubGaussian(α=α, R=cov, n=length(x))
 end
 

test/runtests.jl

@@ -76,20 +76,31 @@ end
     sampletypes = [Float32,Float64]
     stabletypes = [AlphaStable,SymmetricAlphaStable]
     αs = [0.6:0.1:2,1:0.1:2]
+    betas = [-1:0.5:1,0.0]
+    sc = 2.0
     for sampletype ∈ sampletypes
         for (i, stabletype) in enumerate(stabletypes)
             for α in αs[i]
-                d1 = AlphaStable(α=sampletype(α))
-                s = rand(rng,d1, 200000)
-                @test eltype(s) == sampletype
+                for β in betas[i]
+                    d1 = if stabletype == AlphaStable 
+                        stabletype(α=sampletype(α), β=sampletype(β), scale=sampletype(sc))
+                    else
+                        stabletype(α=sampletype(α), scale=sampletype(sc))
+                    end
+                    s = rand(rng, d1, 10^6)
+                    @test eltype(s) == sampletype
 
-                d2 = fit(stabletype, s)
-                @test typeof(d2.α) == sampletype 
+                    d2 = fit(stabletype, s)
+                    @test typeof(d2.α) == sampletype 
 
-                @test d1.α ≈ d2.α rtol=0.1
-                stabletype != SymmetricAlphaStable && @test d1.β ≈ d2.β atol=0.2
-                @test d1.scale ≈ d2.scale rtol=0.1
-                @test d1.location ≈ d2.location atol=0.1
+                    @test d1.α ≈ d2.α rtol=0.1
+                    if (stabletype != SymmetricAlphaStable) && (α != 2)
+                        @test d1.β ≈ d2.β atol=0.2
+                    end
+                    # the quantile method is less accurate
+                    @test d1.scale ≈ d2.scale rtol=0.2 * sc
+                    @test d1.location ≈ d2.location atol=0.9 * sc
+                end
             end
 
             xnormal = rand(rng,Normal(3.0, 4.0), 96000)
@@ -119,7 +130,7 @@ end
         @test d3.α ≈ α rtol=0.2
         @test d3.β ≈ 0 atol=0.2
         @test d3.scale ≈ 1 rtol=0.2
-        @test d3.location ≈ 0 atol=0.1
+        @test d3.location ≈ 0 atol=0.2
     end
 
     d4 = AlphaSubGaussian(α=1.5, n=96000)