# Bug in AlphaStableDistributions?

**URL:** <https://discourse.julialang.org/t/bug-in-alphastabledistributions/99868>\
**Category:** Statistics\
**Tags:** question\
**Created:** [June 5, 2023, 4:46am UTC](https://discourse.julialang.org/t/bug-in-alphastabledistributions/99868 "2023-06-05T04:46:34Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![chrhum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrhum/32/52683_2.png) [@chrhum](https://discourse.julialang.org/u/chrhum)\
**Post date:** [June 5, 2023, 4:46am UTC](https://discourse.julialang.org/t/bug-in-alphastabledistributions/99868/1 "2023-06-05T04:46:34Z")

</div>

I could not reconcile the output of a classical sampling algorithm for the Gumbel Copula when using `AlphaStableDistributions.jl`. I wonder if there is a bug in its simulation algorithm or if I have overseen something?

A classical sampling algorithm for the (nested) Gumbel Copula, see e.g [Hofert et al](https://cran.r-project.org/web/packages/copula/vignettes/nacopula-pkg.pdf) page 4, is based on sampling from a stable distribution. It produces in particular uniformly distributed random numbers `u` in the unit interval:

```plaintext
using AlphaStableDistributions, Statistics

# Select parameters for calculation
theta = 1.2 # needs to be > 1
nsamples = 10^6

# Algorithm 
# Step (1) 
scal(th) = (cos(π / (2.0 * th)))^th
beta = 1.0
s = AlphaStable(α=1.0 / theta, β=beta, scale=scal(theta), location=0.0)

v = rand(s, nsamples)

# Step (2)
r = -log.(rand(nsamples)) 

# Step (3)
psi(x, th) = exp(-x^(1.0 / th))
u = psi.(r ./ v, theta)

```

Assumption: The algorithm is specified for the so-called 1-paramaterization in [Nolan](https://edspace.american.edu/jpnolan/wp-content/uploads/sites/1720/2020/09/Chap1.pdf) which is probably the choice of `AlphaStableDistributions` as it implements the mean, where defined, as the location parameter.

In view of the sample size, `mean` and `var` of `u` are too far away from the theoretical values `0.5` and `0.08333...`

```julia
julia> (mean(u), var(u))
(0.4578483790138697, 0.10688388000233129)

```

Replacing the samples `v` in Step (2) above using the algorithm [Nolan, Theorem 1.19, page 21] the result of the test is as expected:

```plaintext
# sketch for a fix
function stablernd(alpha, beta) # is S(alpha, beta,1,0;1) distributed
    alpha != 1.0 || error("special case missing")
    t0 = atan(beta * tan((pi * alpha) / 2)) / alpha # numerical problems if alpha near 1
    Theta = pi * (rand() - 0.5)
    W = -log(rand())
    term1 = sin(alpha * (t0 + Theta)) / (cos(alpha * t0) * cos(Theta))^(1 / alpha)
    term2 = ((cos(alpha * t0 + (alpha - 1) * Theta)) / W)^((1 - alpha) / alpha)
    return term1 * term2
end

# replace v and recalculate u from above
for k in 1:nsamples
    v[k] = stablernd(1 / theta, beta) * scal(theta)
end
u = psi.(r ./ v, theta)

```

Indeed, …

```julia
julia> (mean(u), var(u))
(0.5003790541200495, 0.08318313400529477)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 5, 2023, 4:49am UTC](https://discourse.julialang.org/t/bug-in-alphastabledistributions/99868/2 "2023-06-05T04:49:26Z")

</div>

Have you opened an issue on the repo? This seems like more of a bug report than a Discourse thread.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [June 5, 2023, 5:32am UTC](https://discourse.julialang.org/t/bug-in-alphastabledistributions/99868/3 "2023-06-05T05:32:56Z")

</div>

This indeed looks like a bug report for AlphaStableDistributions.jl. Note that you could reproduce your first piece of code, including in more dimensions, through :

```julia
using Copulas, Statistics
θ = 1.2
d = 3
G = GumbelCopula(d,θ)
nsamples = 10^6
u = rand(G,nsamples)
(mean(u),var(u))

```

And indeed I get the same weird `(0.4581850334671335, 0.10689564822591738)`, even for dimensions `d > 2`, while I totally agree with you that the mean should be a lot closer to 1/2.

This is using the same AlphaStable parametrisation as you are. Edit: The [`Copulas.jl`’](https://github.com/lrnv/Copulas.jl/tree/main)s one is implemented [in this file](https://github.com/lrnv/Copulas.jl/blob/main/src/univariate_distributions/AlphaStable.jl) so maybe it requires a separate fix from [`AlphaStableDistributions.jl`](https://github.com/org-arl/AlphaStableDistributions.jl).

---

<div class="post-metadata">

**Author:** ![chrhum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrhum/32/52683_2.png) [@chrhum](https://discourse.julialang.org/u/chrhum)\
**Post date:** [June 6, 2023, 6:10am UTC](https://discourse.julialang.org/t/bug-in-alphastabledistributions/99868/4 "2023-06-06T06:10:23Z")

</div>

Thank you for your feedback. I opened an issue on the repo.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [June 6, 2023, 8:51am UTC](https://discourse.julialang.org/t/bug-in-alphastabledistributions/99868/5 "2023-06-06T08:51:58Z")

</div>

So both `Copulas.jl` and `AlphaStableDistributions.jl` were fixed and registered, so that now after updating your packages you have :

```julia
using Copulas, Statistics, Random
Random.seed!(123)
θ = 1.2
d = 2
G = GumbelCopula(d,θ)
nsamples = 10^4
u = rand(G,nsamples)
(mean(u),var(u))

```

returning `(0.5038099870266098, 0.08355682703493884)` as it should be.

Thanks !

PS: May I ask what it is you were trying to do with this ?

---

<div class="post-metadata">

**Author:** ![chrhum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrhum/32/52683_2.png) [@chrhum](https://discourse.julialang.org/u/chrhum)\
**Post date:** [June 7, 2023, 5:01am UTC](https://discourse.julialang.org/t/bug-in-alphastabledistributions/99868/6 "2023-06-07T05:01:24Z")

</div>

Many thanks for this fix and the swift response.

P.S.: I detected the bug when I translated some Matlab code I wrote long ago to Julia in order to sample from nested Gumbel copulas - recursively through the tree with theta parameters in the vertices. I found them quite useful in building actuarial toy models.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [June 7, 2023, 10:00am UTC](https://discourse.julialang.org/t/bug-in-alphastabledistributions/99868/7 "2023-06-07T10:00:09Z")

</div>

> [@chrhum](#):
>
> P.S.: I detected the bug when I translated some Matlab code I wrote long ago to Julia in order to sample from nested Gumbel copulas - recursively through the tree with theta parameters in the vertices. I found them quite useful in building actuarial toy models.

The nested Archimedean trees copulas are on my todo-list for [`Copulas.jl`](https://github.com/lrnv/Copulas.jl). If you have something that works correctly, maybe you could propose it for a PR ? If you dont know how to do such a contribution to an open source repository, you may learn (its a great tool to have) or i could help you.