RNG for alpha stable isn't correct #38

Activity

1. mchitre

   assigned

   tom-plaa

   on Dec 29, 2022
2. tom-plaa commented on Jan 4, 2023

   

   tom-plaa

   on Jan 4, 2023

   Last edited by tom-plaa

   ContributorAuthor

   More actions

   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.
3. tom-plaa commented on Jan 4, 2023

   

   tom-plaa

   on Jan 4, 2023

   Last edited by tom-plaa

   ContributorAuthor

   More actions

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