Assess stable distributed numbers generation

[tail-call]

Is it any good? Does it really work? Are there good implementations that already exist?

See cgtnnlib/NoiseGenerator.py:

def stable_noise_func(alpha, beta, size=1):
    """
    Generates stable noise samples using the Chambers-Mallows-Stuck (CMS) algorithm.

    Args:
      alpha: Stability parameter (0 < alpha <= 2).  alpha=2 corresponds to Gaussian. alpha=1 Cauchy.
      beta: Skewness parameter (-1 <= beta <= 1). beta=0 is symmetric.
      size: Number of samples to generate.

    Returns:
      A numpy array of stable noise samples.

    Raises:
      ValueError: If alpha or beta are outside the allowed ranges.
    """

    if not 0 < alpha <= 2:
        raise ValueError("alpha must be in the range (0, 2]")
    if not -1 <= beta <= 1:
        raise ValueError("beta must be in the range [-1, 1]")

    U = np.random.uniform(-np.pi / 2, np.pi / 2, size=size)
    E = np.random.exponential(1, size=size)

    if alpha != 1:
        term1 = np.sin(alpha * (U + beta * np.pi / 2))
        term2 = np.cos(U) ** (-1 / alpha)
        term3 = np.cos(U - alpha * (U + beta * np.pi / 2))

        # Cast to complex to avoid nans due to exponentiation
        X = term1 * term2 * abs(np.cfloat(E / term3) ** ((1 - alpha) / alpha))

    else:
        X = (2 / np.pi) * (
            ((np.pi / 2) + beta * U) * np.tan(U)
            - beta * np.log((E * np.cos(U)) / ((np.pi / 2) + beta * U))
        )

    return