Merge branch 'master' into compathelper/new_version/2019-12-31-03-05-28-992-167003063

Author: baggepinnen

Project.toml

@@ -16,6 +16,7 @@ ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
 MAT = "0.7"
+StatsBase = "0.32"
 julia = "1.1"
 
 [extras]

README.md

@@ -19,7 +19,7 @@ julia> using AlphaStableDistributions
 julia> d1 = AlphaStable()
 AlphaStable{Float64}(α=1.5, β=0.0, scale=1.0, location=0.0)
 
-julia> s = [rand(d1) for _ in 1:100000];
+julia> s = rand(d1, 100000);
 
 julia> d2 = fit(AlphaStable, s)
 

src/AlphaStableDistributions.jl

@@ -121,26 +121,40 @@ function Base.rand(rng::AbstractRNG, d::AlphaStable)
     bϕ = π/2 + β*ϕ
     x = 2/π * (bϕ*tan(ϕ) - β*log(π/2*w*cosϕ/bϕ))
     α == 1 || (x += β * tan(π*α/2))
-    
+
     return loc + scale*x
 end
 
 
 """
-...
+
+Generate alpha-sub-Gaussian (aSG) random numbers.
+
+The implementation is based on https://github.com/ahmd-mahm/alpha-SGNm/blob/master/asgn.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.)
+
+
 # Arguments
-- `alpha`: characteristic exponent associated with the aSGN(m) process. This is
+- `α`: characteristic exponent associated with the aSGN(m) process. This is
 a scalar input and should lie within `collect(1.1:0.01:1.98)`.
 - `R`: covariance matrix of any adjacent `m+1` samples in an aSGN(m) process.
 The dimension of `R` is equal to `m+1`. It should be a symmetric toeplitz matrix.
 The maximum acceptable size of `R` is `10x10`
 - `n`: number of samples required
-...
+
+# Examples
+```jldoctest
+julia> x = rand(AlphaSubGaussian(n=1000))
+```
 """
 Base.@kwdef struct AlphaSubGaussian{T,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 = 10000
+    n::Int
 end
 
 """
@@ -184,23 +198,10 @@ function subgausscondprobtabulate(α, x1, x2_ind, invRx1, invR, vjoint, nmin, nm
     (1/sqrt(kappa))*kmarg*vjointR/vjointR1
 end
 
-"""
-Generate alpha-sub-Gaussian (aSG) random numbers.
-
-The implementation is based on https://github.com/ahmd-mahm/alpha-SGNm/blob/master/asgn.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.)
-
-# Examples
-```jldoctest
-julia> x = rand(AlphaSubGaussian(n=1000))
-```
-"""
-function Base.rand(rng::AbstractRNG, d::AlphaSubGaussian)
+function Random.rand!(rng::AbstractRNG, d::AlphaSubGaussian, x::AbstractArray)
     α=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`"))
     m = size(R, 1)-1
     funk1 = x -> (2^α)*sin(π*α/2)*gamma((α+2)/2)*gamma((α+x)/2)/(gamma(x/2)*π*α/2)
@@ -226,13 +227,13 @@ function Base.rand(rng::AbstractRNG, d::AlphaSubGaussian)
     S = S/sqrt(sum(abs2,S))
     xtmp = ((sigrootx1*sqrt(A*T))*S)'
     if n<=m
-        x = xtmp[1:n]
+        copyto!(x, @view(xtmp[1:n]))
     else
-        x = zeros(n)
+        # x = zeros(n)
         x[onetom] = xtmp
         vstud = α+m
         norms = pdf(TDist(vstud), 0.0)
-        for i = m+1:n
+        @inbounds for i = m+1:n
             x1 = SVector{m,Float64}(view(x,i-m:i-1))
             mode = modefactor*x1
             norm1 = subgausscondprobtabulate(α, x1, mode, invRx1, invR, vjoint, nmin, nmax, step, rind, kappa, k1, k2, kmarg)
@@ -252,4 +253,7 @@ function Base.rand(rng::AbstractRNG, d::AlphaSubGaussian)
     x
 end
 
+
+Base.rand(rng::AbstractRNG, d::AlphaSubGaussian) = rand!(rng, d, zeros(d.n))
+
 end # module

test/runtests.jl

@@ -4,7 +4,8 @@ using Test
 @testset "AlphaStableDistributions.jl" begin
 
 d1 = AlphaStable()
-s = [rand(d1) for _ in 1:100000]
+s = rand(d1, 100000)
+
 d2 = fit(AlphaStable, s)
 
 @test d1.α ≈ d2.α rtol=0.1
@@ -13,8 +14,17 @@ d2 = fit(AlphaStable, s)
 @test d1.location ≈ d2.location atol=0.03
 
 
-x = rand(AlphaSubGaussian(n=96000))
+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
 
 end
 # 362.499 ms (4620903 allocations: 227.64 MiB)