1"""Gaussian compensation for the small jumps of a stable Levy noise."""
 2import math
 3import numpy as np
 4from scipy.stats import levy_stable
 5
 6
 7def small_jump_variance(eps, alpha, c=1.0):
 8    return c * eps ** (2.0 - alpha) / (2.0 - alpha)
 9
10
11def small_jump_third_abs_moment(eps, alpha, c=1.0):
12    return c * eps ** (3.0 - alpha) / (3.0 - alpha)
13
14
15def stable_scale(alpha, c=1.0):
16    # Symmetric Levy measure with total radial intensity c*u^(-1-alpha)du.
17    # Integral (cos(tx)-1)c*x^(-1-alpha)dx = -scale**alpha*|t|**alpha.
18    if abs(alpha - 1.0) < 1e-12:
19        return c * math.pi / 2.0
20    return -c * math.gamma(-alpha) * math.cos(math.pi * alpha / 2.0)
21
22
23def sample_exact_small_symmetric(eps, alpha, n, rng, c=1.0):
24    """Sample the centered residual over 0<u<eps.
25
26    It is represented as a symmetric alpha-stable variable minus its jumps
27    above eps. This avoids imposing an artificial lower cutoff.
28    """
29    scale = stable_scale(alpha, c) ** (1.0 / alpha)
30    # scipy's random_state accepts a RandomState, not Generator.
31    rs = np.random.RandomState(rng.integers(0, 2**31 - 1))
32    full = levy_stable.rvs(alpha, 0.0, loc=0.0, scale=scale,
33                           size=n, random_state=rs)
34    rate = c * eps ** (-alpha) / alpha
35    k = rng.poisson(rate, size=n)
36    total = int(k.sum())
37    if total:
38        # Conditional radial law: P(U>u)=(u/eps)^(-alpha).
39        u = eps * rng.random(total) ** (-1.0 / alpha)
40        signs = rng.choice(np.array([-1.0, 1.0]), size=total)
41        owners = np.repeat(np.arange(n), k)
42        large = np.bincount(owners, weights=signs * u, minlength=n)
43    else:
44        large = np.zeros(n)
45    return full - large, int(total)
46
47
48def approximate_samples(eps, alpha, n, rng, c=1.0):
49    exact, jumps = sample_exact_small_symmetric(eps, alpha, n, rng, c)
50    naive = np.zeros(n)
51    gaussian = rng.normal(0.0, math.sqrt(small_jump_variance(eps, alpha, c)), n)
52    return exact, naive, gaussian, jumps
53
54
55def w1_1d(x, y):
56    """Empirical one-dimensional W1 for equal-size samples."""
57    return float(np.mean(np.abs(np.sort(x) - np.sort(y))))