Merge branch 'master' of github.com:org-arl/AlphaStableDistributions.jl

Author: baggepinnen

.travis.yml

@@ -2,7 +2,6 @@
 language: julia
 os:
   - linux
-  - osx
 julia:
   - 1.3
   - nightly

Project.toml

@@ -15,13 +15,13 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
-Distributions = "0.21"
+Distributions = "0.21, 0.22"
 MAT = "0.7"
-StatsBase = "0.32"
-julia = "1.1"
-SpecialFunctions = "0.8, 0.9"
+SpecialFunctions = "0.8, 0.9, 0.10"
 StaticArrays = "0.12"
+StatsBase = "0.32"
 ToeplitzMatrices = "0.6"
+julia = "1.1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -5,7 +5,7 @@
 
 This library is a port of functionality from [arlpy](https://github.com/org-arl/arlpy/blob/master/arlpy/stable.py). The two distributions supported are
 - [alpha-stable distribution](https://en.wikipedia.org/wiki/Stable_distribution) (`rand` and `fit`)
-- [alpha sub-Gaussian distribution with memory](https://arl.nus.edu.sg/twiki6/pub/ARL/BibEntries/SigProc2016RandomVariate.pdf) (`rand`)
+- [alpha sub-Gaussian distribution with memory](https://arl.nus.edu.sg/twiki6/pub/ARL/BibEntries/SigProc2016RandomVariate.pdf) (`rand` and `fit`)
 
 ## Installation
 ```julia

src/AlphaStableDistributions.jl

@@ -76,10 +76,10 @@ function Distributions.fit(d::Type{<:AlphaStable}, x)
     else
         α = 0.
         j = findfirst(>=(an), _ena) # _np.where(an <= _ena[:,0])[0]
-        if j !== nothing
-            t = (an-_ena[j,1])/(_ena[j+1]-_ena[j])
-            α = (21-j-t)/10
-        end
+        (j === nothing || j == length(_ena)) && (j = length(_ena)-1)
+        t = (an-_ena[j])/(_ena[j+1]-_ena[j])
+        α = (21-j-t)/10
+
     end
     if α < 0.5
         α = 0.5
@@ -185,15 +185,15 @@ function subgausscondprobtabulate(α, x1, x2_ind, invRx1, invR, vjoint, nmin, nm
     if r<nmin
         grad = (vjoint[m+1, 1]-k2[2])/nmin
         cons = k2[2]
-        vjointR = grad*r1+cons
+        vjointR = grad*r+cons
     elseif r>nmax
-        vjointR = α*k1[2]*(r^(-α-m))
+        vjointR = α*k1[2]*(r^(-α-m-1))
     else
         ti = (log10(r)-log10(nmin))/step+1
         tempind = (floor(Int, ti), ceil(Int, ti))
         grad = (vjoint[m+1, tempind[1]]-vjoint[m+1, tempind[2]])/(rind[tempind[1]]-rind[tempind[2]])
         cons = vjoint[m+1, tempind[1]]-grad*rind[tempind[1]]
-        vjointR = grad*r1+cons
+        vjointR = grad*r+cons
     end
     (1/sqrt(kappa))*kmarg*vjointR/vjointR1
 end
@@ -256,4 +256,41 @@ end
 
 Base.rand(rng::AbstractRNG, d::AlphaSubGaussian) = rand!(rng, d, zeros(d.n))
 
+"""
+    fit(d::Type{<:AlphaSubGaussian}, x, m; p=1.0)
+
+Fit an aSGN(m) model to data via the covariation method.
+
+The covariation method requires an additional parameter `p`. Ideally, 1 < p < α. In most practical impulsive scenarios p=1.0 is sufficient.
+`m` is the number of lags in the covariance matrix.
+
+The implementation is based on https://github.com/ahmd-mahm/alpha-SGNm/blob/master/param_est/asgnfit.m
+
+Reference:
+A. Mahmood and M. Chitre, "Generating random variates for stable sub-Gaussian processes
+with memory", Signal Processing, Volume 131, Pages 271-279, 2017.
+(https://doi.org/10.1016/j.sigpro.2016.08.016.)
+"""
+function Distributions.fit(d::Type{<:AlphaSubGaussian}, x::AbstractVector{T}, m::Integer; p=one(T)) where T
+    d1   = fit(AlphaStable, x)
+    α    = d1.α; scale=d1.scale
+    cov  = zeros(T, m+1, m+1)
+    xlen = length(x)
+    c    = ((sum(abs.(x).^p)/xlen)^(1/p))/scale
+    for i in 1:m
+        tempxlen = xlen-mod(xlen, i)
+        xtemp = reshape(x[1:end-mod(xlen, i)], i, tempxlen÷i)
+        if mod(tempxlen÷i, 2) != 0
+            xtemp = xtemp[:, 1:end-1]
+            tempxlen = size(xtemp, 1)*size(xtemp, 2)
+        end
+        xtemp = reshape(xtemp', 2, tempxlen÷2)
+        r = (2/(c^p))*(scale^(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
+    AlphaSubGaussian(α=α, R=cov, n=length(x))
+end
+
 end # module

test/runtests.jl

@@ -3,28 +3,39 @@ using Test, Random
 
 @testset "AlphaStableDistributions.jl" begin
 
-d1 = AlphaStable()
-s = rand(d1, 100000)
+    for α in 0.5:0.1:2
+        d1 = AlphaStable(α=α)
+        s = rand(d1, 100000)
 
-d2 = fit(AlphaStable, s)
+        d2 = fit(AlphaStable, s)
 
-@test d1.α ≈ d2.α rtol=0.1
-@test d1.β ≈ d2.β rtol=0.1
-@test d1.scale ≈ d2.scale rtol=0.1
-@test d1.location ≈ d2.location atol=0.03
+        @test d1.α ≈ d2.α rtol=0.1
+        @test d1.β ≈ d2.β rtol=0.1
+        @test d1.scale ≈ d2.scale rtol=0.1
+        @test d1.location ≈ d2.location atol=0.03
+    end
 
+    for α in 1.1:0.1:1.9
 
-d = AlphaSubGaussian(n=96000)
-x = rand(d)
-x2 = copy(x)
-rand!(d, x2)
-@test x != x2
+        d = AlphaSubGaussian(n=96000, α=α)
+        x = rand(d)
+        x2 = copy(x)
+        rand!(d, x2)
+        @test x != x2
 
-d3 = fit(AlphaStable, x)
-@test d3.α ≈ 1.5 rtol=0.2
-@test d3.β == 0
-@test d3.scale ≈ 1 rtol=0.2
-@test d3.location ≈ 0 atol=0.03
+        d3 = fit(AlphaStable, x)
+        @test d3.α ≈ α rtol=0.2
+        @test d3.β == 0
+        @test d3.scale ≈ 1 rtol=0.2
+        @test d3.location ≈ 0 atol=0.03
+    end
+
+    d4 = AlphaSubGaussian(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)