Correct alphastable rng for alpha in (1,2)

[tom-plaa]

Trying to fix #38.

[ymtoo]

@tom-plaa rand(rng::AbstractRNG, d::AlphaStable{T}) is currently based on McCulloch's implementation (1996). Can you share the reference for the proposed modification?

[tom-plaa]

(I replied by email but I don't think it worked, so copying it here)

The references I based myself on are the original paper by Chambers, Mallows and Stuck (https://doi.org/10.2307/2285309) and the survey by Kawai and Masuda (https://doi.org/10.1016/j.cam.2010.12.014).

For a direct reference, see for example formula (2.2) on Kawai-Masuda and remember that the parametrization used in that paper is the one with the Lévy measure, so we have (using their notation):

scale^alpha = - a * gamma(-alpha) * cos(pi*alpha/2)

Additionally, see formula (2.3) on Chamber-Mallows-Stuck (the original reference for the algorithm).

NOTE: To arrive at the formula that is equivalent to the one here in the code from Kawai-Masuda's, we need to make a few simple manipulations using the properties of sin(a+b), sin(a-b), cos(a+b) and cos(a-b). I can post them here if needed.

[ymtoo]

Much appreciated if you can show briefly the manipulation of how to arrive at the formula.

And can you add test cases (runtests.jl) such that the current implementation throws an error and the modified code passes the test?

[tom-plaa]

Regarding the formula, I will provide the derivation in the next following days; I will scan my sheets as it would be very tiresome to type that all out in LaTeX.

About the tests, it seems the current ones are based on the fit procedure and both schemes pass successfully. The test case I presented in the issue clearly shows us different behaviors for the two schemes, where this scheme seems to follow the numerically inverted pdf.
The question of how to test so that one passes and the other fails seems to imply we need to make use of the numerical invertion of the characteristic function to get back the cdf and then compare some distribution distance (kolmogorov-smirnov, anderson-darling...) using the numerical cdf and montecarlo generated samples, and then one scheme should have a noticeably smaller distance. The problem here is that this implies a new dependency, the module I used for this in #38 was my own code that isn't in the Julia registry (https://gitlab.com/tom.plaa/characteristicinvfourier.jl).
How should we include this functionality? We could include a new function in the test script that does the invertion for a given object of type AlphaStable, using the cf defined in the module, but it will inevitably have to include all the steps required for numerical invertion of the cf and the usage of the FFTW package.

[ymtoo]

@tom-plaa take your time.

CharacteristicInvFourier.jl looks useful. Do you plan to register the package? If yes, we can add the package as test-dependencies. One of the potential test cases could be as follows.

julia> using HypothesisTests, Test

julia> test_pkg_cdf = cdf_invfourier_fft(test_pkg_cf, npower=18, xbounds=(-10.0, 1000));

julia> Distributions.cdf(::AlphaStable, x) = test_pkg_cdf(x)

julia> t1 = ExactOneSampleKSTest(test_pkg_mcsamples, test_stable)
Exact one sample Kolmogorov-Smirnov test
----------------------------------------
Population details:
    parameter of interest:   Supremum of CDF differences
    value under h_0:         0.0
    point estimate:          0.296395

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           <1e-99

Details:
    number of observations:   100000000

julia> t2 = ExactOneSampleKSTest(test_new_mcsamples, test_stable)
Exact one sample Kolmogorov-Smirnov test
----------------------------------------
Population details:
    parameter of interest:   Supremum of CDF differences
    value under h_0:         0.0
    point estimate:          0.00012986

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.0686

Details:
    number of observations:   100000000

julia> @test t1.δ < 1e-3 # supremum of CDF differences
Test Failed at REPL[52]:1
  Expression: t1.δ < 0.001
   Evaluated: 0.296394601850442 < 0.001
ERROR: There was an error during testing

julia> @test t2.δ < 1e-3
Test Passed

[tom-plaa]

I have nothing against registering CharacteristicInvFourier.jl, it's just that I fear it isn't robust enough to be a registered package because it has no testing pipeline; and unfortunately I do not have the time to develop one at the moment. I will probably ask for some opinions back at the Julia discourse to see if it's generally considered acceptable to register it as it is.

Regarding the tests you proposed, that's exactly what I had in mind.

[tom-plaa]

Hello, sorry I haven't yet updated this, I am currently busy with other topics. I will try to find an opportunity at least in February (if not sooner).

[tom-plaa]

I have gotten the CharacteristicInvFourier package automerged on the julia registry. If it gets approved I will do the missing things in this PR.

[tom-plaa]

Here are my computations. I will now work on the unit tests.

[tom-plaa]

After analyzing this further, I am not sure that the formula I provided is correct for general beta (it is for beta=1). There's the paper with Weron's algorithm at
On the Chambers-Mallows-Stuck method for simulating skewed stable random variables which I will be reading over to make sure we have the correct formula for general beta.