RNG for alpha stable isn't correct

[tom-plaa]

I was doing some numerical tests where I ended up comparing rng results from this package by using the characteristic function of the target distribution. At least for alpha in (1,2), the results here do not coincide with the results from inverting the cdf.

It has to do with the (1-alpha)/alpha exponent not being applied to the trignometric functions and also maybe with the generation of the exponential random variable. The formula here is definitely not correct.

I have made a few tests on my machine with some corrections and they seem to work, so I will submit a fix in the near future and we can evaluate the results.

I created this issue to raise awareness, as I plan to submit a fix to this in the next few days/weeks.

[tom-plaa]

To expand a bit on the issue description, here is a script that elucidates what I meant. The plot is truncated so we can see the curves near the mode better, but the rng generates gigantic samples a lot more frequently than theoretically. Calculating any distribution distance metric (kolmogorov-smirnov, anderson-darling) will yield enormous distances as well.

using Distributions, AlphaStableDistributions
using SpecialFunctions
using Plots, StatsPlots
using CharacteristicInvFourier # https://gitlab.com/tom.plaa/characteristicinvfourier.jl

test_stable = AlphaStable(1.5, 1.0, 1.0, 0.0)

test_pkg_mcsamples = rand(test_stable, 10^8)

test_pkg_cf(x) = cf(test_stable, x)

test_pkg_pdf = pdf_invfourier_fft(test_pkg_cf, npower=18, xbounds=(-10.0, 1000))

xrange = -10:0.01:100

plot(xrange, test_pkg_pdf.(xrange))
histogram!(test_pkg_mcsamples, normalize=:pdf, bins=10^7)
plot!(xlims=(-2.0, 10.0))

Afterwards, I will try to generate samples with what I think is the corrected formula, and post the same type of comparison as well.

[tom-plaa]

Here is the proposed version, using the same test case as above and the modified formula:

using Distributions, AlphaStableDistributions
using SpecialFunctions
using Plots, StatsPlots
using CharacteristicInvFourier # https://gitlab.com/tom.plaa/characteristicinvfourier.jl
using Random

test_stable = AlphaStable(1.5, 1.0, 1.0, 0.0)

# test_pkg_mcsamples = rand(test_stable, 10^8)

test_pkg_cf(x) = cf(test_stable, x)

test_pkg_pdf = pdf_invfourier_fft(test_pkg_cf, npower=18, xbounds=(-10.0, 1000))

# only the case α ∈ (1.0, 2.0)
function new_rng(rng::AbstractRNG, d::AlphaStable{T}) where {T<:AbstractFloat}
    α=d.α; β=d.β; scale=d.scale; loc=d.location

    ϕ = (rand(rng, T) - 0.5) * π

    w = randexp(rng, T)

    # β == zero(T) && (return loc + scale * ((cos((1-α)*ϕ) / w)^(one(T)/α - one(T)) * sin(α * ϕ) / cos(ϕ)^(one(T)/α)))

    cosϕ = cos(ϕ)
    ζ = β * tan(π * α / 2)
    aϕ = α * ϕ
    a1ϕ = (one(T) - α) * ϕ

    unit_term = ((sin(aϕ) + ζ*cos(aϕ)) / (cosϕ^(1/α))) * ((cos(a1ϕ) + ζ*sin(a1ϕ)) / w) ^ ((1-α)/α)

    return loc + scale * unit_term
end

test_rng=Xoshiro()

test_new_mcsamples = map(x->new_rng(test_rng, test_stable), zeros(10^8))

xrange = -10:0.01:100

plot(xrange, test_pkg_pdf.(xrange), label="numerical pdf")
# histogram!(test_pkg_mcsamples, normalize=:pdf, bins=10^7, label="pkg montecarlo")
histogram!(test_new_mcsamples, normalize=:pdf, bins=10^7, label="new montecarlo")
plot!(xlims=(-2.0, 10.0))