Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

RNG for alpha stable isn't correct #38

Open
tom-plaa opened this issue Dec 28, 2022 · 2 comments · May be fixed by #39
Open

RNG for alpha stable isn't correct #38

tom-plaa opened this issue Dec 28, 2022 · 2 comments · May be fixed by #39
Assignees
@tom-plaa

Comments

@tom-plaa
Copy link
Contributor

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
Copy link
Contributor Author

tom-plaa commented Jan 4, 2023 •
edited
Loading

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))

numericalpdf_vs_montecarlohistogram

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
Copy link
Contributor Author

tom-plaa commented Jan 4, 2023 •
edited
Loading

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))

numericalpdf_vs_newsamples

@tom-plaa tom-plaa linked a pull request Jan 4, 2023 that will close this issue
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@tom-plaa tom-plaa

Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Correct alphastable rng for alpha in (1,2) tom-plaa/AlphaStableDistributions.jl
1 participant
@tom-plaa