1"""Exact exponential-model identities (Theorem 3 ingredients) checked numerically on random inputs."""2import numpy as np3from scipy import stats4from pamle.jackknife import vandermonde_weights, leading_bias_coefficient567def test_gamma_inverse_moments():8 rng = np.random.default_rng(1)9 for n in (5, 17, 64):10 theta0 = 2.511 S = rng.gamma(n, 1 / theta0, size=400000)12 th = n / S13 assert abs(th.mean() - theta0 * n / (n - 1)) < 0.01 * theta014 assert abs(th.var() - theta0 ** 2 * n ** 2 / ((n - 1) ** 2 * (n - 2))) < 0.05 * theta0 ** 2 * n ** 2 / ((n - 1) ** 2 * (n - 2))151617def test_newton_identity():18 rng = np.random.default_rng(2)19 for _ in range(1000):20 x = rng.uniform(0.1, 5); xbar = rng.uniform(0.1, 5); th = 1 / xbar21 g = 1 / x - xbar; H = -1 / x ** 222 newton = x - g / H23 assert abs(newton - (2 * x - x * x * xbar)) < 1e-1224 assert abs((newton - th) + (x - th) ** 2 / th) < 1e-12252627def test_newton_basin():28 th = 1.329 for x in (0.01, 0.5, 1.0, 2.0, 2.59):30 e = x - th31 for _ in range(60):32 e = -e * e / th33 assert abs(e) < 1e-634 for x in (-0.1, 2.61, 3.0):35 e = x - th36 for _ in range(12):37 e = -e * e / th38 assert abs(e) > 139 e = 0.0 + (2 * th - th) # x = 2 theta_hat -> maps to 0 -> stays at 0 (error -theta_hat)40 e = -e * e / th41 assert abs(e + th) < 1e-12424344def test_exact_bias_of_split_estimators():45 """E[J_r] - theta0 = theta0 * sum_j w_j s_j/(m - s_j) exactly, leading term beta_r theta0 / m^r."""46 theta0 = 1.047 for s in [(1, 2), (1, 2, 4), (1, 2, 4, 8)]:48 w = vandermonde_weights(s); r = len(s)49 for m in (64, 256, 1024):50 bias = sum(wj * theta0 * m / (m - sj) for wj, sj in zip(w, s)) - theta051 if s == (1, 2):52 assert abs(bias + 2 * theta0 / ((m - 1) * (m - 2))) < 1e-1453 if s == (1, 2, 4):54 assert abs(bias - 8 * theta0 / ((m - 1) * (m - 2) * (m - 4))) < 1e-1455 assert abs(bias * m ** r / (float(leading_bias_coefficient(s)) * theta0) - 1) < 20 / m565758def test_boundary_constant_numeric():59 """Deterministic version of Theorem 3: plug e0 = beta_r theta0 c^r N^{-1/2^(q+1)} into the q-fold Newton identity."""60 theta0, c = 1.7, 0.861 for s, q in [((1, 2), 1), ((1, 2, 4), 1), ((1, 2), 2), ((1, 2, 4, 8), 1)]:62 r = len(s); beta = float(leading_bias_coefficient(s))63 for N in (1e12, 1e16):64 alpha = 1 - 1 / (2 ** (q + 1) * r); m = N ** (1 - alpha) / c65 e = beta * theta0 * m ** (-r)66 for _ in range(q):67 e = -e * e / theta068 lim = -theta0 * (beta * c ** r) ** (2 ** q)69 assert abs(np.sqrt(N) * e - lim) < 1e-8 * abs(lim)70 assert abs((-theta0 * ((-2) * c ** 2) ** 2) - (-4 * theta0 * c ** 4)) < 1e-14717273def test_exact_variance_of_half_block_jackknife():74 """Var(J2) = theta0^2 m^2 (m^2 - 6m + 17)/((m-4)(m-2)^2(m-1)^2) (SymPy, theory/symbolic_checks.py S6) vs simulation."""75 rng = np.random.default_rng(5); m = 32; B = 40000076 halves = rng.gamma(m // 2, 1.0, size=(B, 2))77 full = m / halves.sum(axis=1); half = (m / 2) / halves78 J2 = 2 * full - half.mean(axis=1)79 exact = m ** 2 * (m ** 2 - 6 * m + 17) / ((m - 4) * (m - 2) ** 2 * (m - 1) ** 2)80 assert abs(J2.var() / exact - 1) < 0.0281 assert abs(J2.mean() - (1 - 2 / ((m - 1) * (m - 2)))) < 0.00282