Infinite random variates

[mchitre]

rand(SymmetricAlphaStable(1.75, 1f0, 0f0)) sometimes generates infinite return value.

[ymtoo]

I can reproduce the error using

julia> for _ ∈ 1:1000000000
         x = rand(SymmetricAlphaStable(1.75f0, 1f0, 0f0))
         isinf(x) && throw("Infinity")
       end

ERROR: "Infinity"
Stacktrace:
 [1] top-level scope
   @ ./REPL[19]:3

[ymtoo]

To simulate w, which is an independent exponential random variable with mean 1 (https://en.wikipedia.org/wiki/Stable_distribution#Simulation_of_stable_variates), x=rand(rng, T) is used to produce $x\in[0,1)$. This is not correct as $x$ should be uniform on $(0, 1)$.

When x=0, w=-log(0)=Inf and (cos((one(T)-α)*ϕ) / w)^(one(T)/α - one(T)) yields Inf since (one(T)/α - one(T)) is a negative value.