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 page 4, is based on sampling from a stable distribution. It produces in particular uniformly distributed random numbers `u`

in the unit interval:

```
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 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> (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:

```
# 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> (mean(u), var(u))
(0.5003790541200495, 0.08318313400529477)
```