Merge pull request org-arl#33 from org-arl/patch-1

fix: rand! for AlphaSubGaussian

Author: ymtoo

src/AlphaStableDistributions.jl

@@ -14,6 +14,8 @@ Base.@kwdef struct AlphaStable{T} <: Distributions.ContinuousUnivariateDistribut
     location::T = zero(α)
 end
 
+AlphaStable(α::Integer, β::Integer, scale::Integer, location::Integer) = AlphaStable(float(α), float(β), float(scale), float(location))
+
 
 # sampler(d::AlphaStable) = error("Not implemented")
 # pdf(d::AlphaStable, x::Real) = error("Not implemented")
@@ -186,7 +188,7 @@ skewness parameter, scale parameter (dispersion^1/α) and location parameter res
 
 α, β, c and δ are computed based on McCulloch (1986) fractile.
 """
-function Distributions.fit(::Type{<:AlphaStable}, x, alg=QuickSort)
+function Distributions.fit(::Type{<:AlphaStable}, x::AbstractArray{T}, alg=QuickSort) where {T<:AbstractFloat}
     sx = sort(x, alg=alg)
     p = quantile.(Ref(sx), (0.05, 0.25, 0.28, 0.5, 0.72, 0.75, 0.95), sorted=true)
     να = (p[7]-p[1]) / (p[6]-p[2])
@@ -210,7 +212,7 @@ function Distributions.fit(::Type{<:AlphaStable}, x, alg=QuickSort)
     else
         δ = ζ - β * c * tan(π*α/2)
     end
-    return AlphaStable(α=α, β=β, scale=c, location=oftype(α, δ))
+    return AlphaStable(α=T(α), β=T(β), scale=T(c), location=T(δ))
 end
 
 Base.@kwdef struct SymmetricAlphaStable{T} <: Distributions.ContinuousUnivariateDistribution
@@ -233,12 +235,12 @@ returns `SymmetricAlphaStable`
 scale is computed based on Fama & Roll (1971) fractile.
 location is the 50% trimmed mean of the sample.
 """
-function Distributions.fit(::Type{<:SymmetricAlphaStable}, x, alg=QuickSort)
+function Distributions.fit(::Type{<:SymmetricAlphaStable}, x::AbstractArray{T}, alg=QuickSort) where {T<:AbstractFloat}
     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
-    an = (p[6]-p[1])/(p[5]-p[2])
+    c = (p[4]-p[3]) / 1.654
+    an = (p[6]-p[1]) / (p[5]-p[2])
     if an < 2.4388
         α = 2.
     else
@@ -251,7 +253,7 @@ function Distributions.fit(::Type{<:SymmetricAlphaStable}, x, alg=QuickSort)
     if α < 0.5
         α = 0.5
     end
-    return SymmetricAlphaStable(α=α, scale=c, location=oftype(α, δ))
+    return SymmetricAlphaStable(α=T(α), scale=T(c), location=T(δ))
 end
 
 """
@@ -267,31 +269,33 @@ 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)
+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) - 0.5) * π
-    if α == 1 && β == 0
-        return loc + scale*tan(ϕ)
+    ϕ = (rand(rng, T) - 0.5) * π
+    if α == one(T) && β == zero(T)
+        return loc + scale * tan(ϕ)
     end
-    w = -log(rand(rng))
+    w = -log(rand(rng, T))
     α == 2 && (return loc + 2*scale*sqrt(w)*sin(ϕ))
-    β == 0 && (return loc + scale * ((cos((1-α)*ϕ) / w)^(1.0/α - 1) * sin(α * ϕ) / cos(ϕ)^(1.0/α)))
+    β == zero(T) && (return loc + scale * ((cos((1-α)*ϕ) / w)^(one(T)/α - one(T)) * sin(α * ϕ) / cos(ϕ)^(one(T)/α)))
     cosϕ = cos(ϕ)
-    if abs(α-1) > 1e-8
-        ζ = β * tan(π*α/2)
+    if abs(α - one(T)) > 1e-8
+        ζ = β * tan(π * α / 2)
         aϕ = α * ϕ
-        a1ϕ = (1-α) * ϕ
-        return loc + scale * (( (sin(aϕ)+ζ*cos(aϕ))/cosϕ * ((cos(a1ϕ)+ζ*sin(a1ϕ))) / ((w*cosϕ)^((1-α)/α)) ))
+        a1ϕ = (one(T) - α) * ϕ
+        return loc + scale * (( (sin(aϕ) + ζ * cos(aϕ))/cosϕ * ((cos(a1ϕ) + ζ*sin(a1ϕ))) / ((w*cosϕ)^((1-α)/α)) ))
     end
     bϕ = π/2 + β*ϕ
-    x = 2/π * (bϕ*tan(ϕ) - β*log(π/2*w*cosϕ/bϕ))
-    α == 1 || (x += β * tan(π*α/2))
+    x = 2/π * (bϕ * tan(ϕ) - β * log(π/2*w*cosϕ/bϕ))
+    α == one(T) || (x += β * tan(π*α/2))
 
-    return loc + scale*x
+    return loc + scale * x
 end
 
+Base.eltype(::Type{<:AlphaStable{T}}) where {T<:AbstractFloat} = T
+
 
 """
 
@@ -318,7 +322,7 @@ The maximum acceptable size of `R` is `10x10`
 julia> x = rand(AlphaSubGaussian(n=1000))
 ```
 """
-Base.@kwdef struct AlphaSubGaussian{T,M<:AbstractMatrix} <: Distributions.ContinuousUnivariateDistribution
+Base.@kwdef struct AlphaSubGaussian{T<:AbstractFloat,M<:AbstractMatrix} <: Distributions.ContinuousUnivariateDistribution
     α::T = 1.50
     R::M = SMatrix{5,5}(collect(SymmetricToeplitz([1.0000, 0.5804, 0.2140, 0.1444, -0.0135])))
     n::Int
@@ -366,7 +370,7 @@ function subgausscondprobtabulate(α, x1, x2_ind, invRx1, invR, vjoint, nmin, nm
 end
 
 
-function Random.rand!(rng::AbstractRNG, d::AlphaSubGaussian, x::AbstractArray)
+function Random.rand!(rng::AbstractRNG, d::AlphaSubGaussian{T}, x::AbstractArray{T}) where {T<:AbstractFloat}
     α=d.α; R=d.R; n=d.n
     length(x) >= n || throw(ArgumentError("length of x must be at least n"))
     α ∈ 1.10:0.01:1.98 || throw(DomainError(α, "α must lie within `1.10:0.01:1.98`"))
@@ -388,11 +392,11 @@ function Random.rand!(rng::AbstractRNG, d::AlphaSubGaussian, x::AbstractArray)
     nmax, nmin, res, rind, vjoint = matdict["Nmax"]::Float64, matdict["Nmin"]::Float64, matdict["res"]::Float64, vec(matdict["rind"])::Vector{Float64}, matdict["vJoint"]::Matrix{Float64}
     step = (log10(nmax)-log10(nmin))/res
     m>size(vjoint, 1)-1 && throw(DomainError(R, "The dimensions of `R` exceed the maximum possible 10x10"))
-    A = rand(AlphaStable(α/2, 1.0, 2*cos(π*α/4)^(2.0/α), 0.0))
-    T = rand(Chisq(m))
+    A = rand(AlphaStable(T(α/2), one(T), T(2*cos(π*α/4)^(2.0/α)), zero(T)))
+    CT = rand(Chisq(m))
     S = randn(m)
     S = S/sqrt(sum(abs2,S))
-    xtmp = ((sigrootx1*sqrt(A*T))*S)'
+    xtmp = ((sigrootx1*sqrt(A*CT))*S)'
     if n<=m
         copyto!(x, @view(xtmp[1:n]))
     else
@@ -421,7 +425,8 @@ function Random.rand!(rng::AbstractRNG, d::AlphaSubGaussian, x::AbstractArray)
 end
 
 
-Base.rand(rng::AbstractRNG, d::AlphaSubGaussian) = rand!(rng, d, zeros(d.n))
+Base.rand(rng::AbstractRNG, d::AlphaSubGaussian) = rand!(rng, d, zeros(eltype(d), d.n))
+Base.eltype(::Type{<:AlphaSubGaussian}) = Float64
 
 """
     fit(d::Type{<:AlphaSubGaussian}, x, m; p=1.0)

test/runtests.jl

@@ -3,38 +3,43 @@ using Test, Random, Distributions
 
 @testset "AlphaStableDistributions.jl" begin
 
+    sampletypes = [Float32,Float64]
     stabletypes = [AlphaStable,SymmetricAlphaStable]
     αs = [0.6:0.1:2,1:0.1:2]
-    for (i, stabletype) in enumerate(stabletypes)
-        for α in αs[i]
-            d1 = AlphaStable(α=α)
-            s = rand(d1, 100000)
+    for sampletype ∈ sampletypes
+        for (i, stabletype) in enumerate(stabletypes)
+            for α in αs[i]
+                d1 = AlphaStable(α=sampletype(α))
+                s = rand(d1, 200000)
+                @test eltype(s) == sampletype
 
-            d2 = fit(stabletype, s)
+                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
-        end
+                @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
+            end
 
-        xnormal = rand(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
+            xnormal = rand(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)
-        d = fit(stabletype, xcauchy)
-        @test d.α ≈ 1 rtol=0.2
-        stabletype != SymmetricAlphaStable && @test d.β ≈ 0 atol=0.2
-        @test d.scale ≈ 4 rtol=0.2
-        @test d.location ≈ 3 rtol=0.1
+            xcauchy = rand(Cauchy(3.0, 4.0), 96000)
+            d = fit(stabletype, xcauchy)
+            @test d.α ≈ 1 rtol=0.2
+            stabletype != SymmetricAlphaStable && @test d.β ≈ 0 atol=0.2
+            @test d.scale ≈ 4 rtol=0.2
+            @test d.location ≈ 3 rtol=0.1
+        end
     end
 
     for α in 1.1:0.1:1.9
-        d = AlphaSubGaussian(n=96000, α=α)
+        d = AlphaSubGaussian(α=α, n=96000)
         x = rand(d)
         x2 = copy(x)
         rand!(d, x2)
@@ -47,14 +52,15 @@ using Test, Random, Distributions
         @test d3.location ≈ 0 atol=0.1
     end
 
-    d4 = AlphaSubGaussian(n=96000)
-    m = size(d4.R, 1)-1
+    d4 = AlphaSubGaussian(α=1.5, n=96000)
+    m = size(d4.R, 1) - 1
     x = rand(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