1#!/usr/bin/env python32"""Symbolic verification of the algebraic identities used in WP28 (SymPy).34Every check raises AssertionError on failure. Run: python3 theory/symbolic_checks.py5Items (numbering follows theory/THEOREM_AUDIT.md):6 S1 Vandermonde extrapolation weights for several split patterns; moment constraints.7 S2 Exponential(rate) model: Gamma inverse moments, E[1/Xbar_n], Var, remainder R_n moments.8 S3 Exact bias of the split estimators J_r in the exponential model and its leading coefficient beta_r.9 S4 Exact Newton identity N(x) = 2x - x^2 Xbar and the error identity N(x) - theta_hat = -(x - theta_hat)^2 / theta_hat.10 S5 q-fold Newton error identity and the boundary limit -theta0 (beta_r c^r)^(2^q).11 S6 Exact variance of the half-block jackknife J_2 in the exponential model (nested dependence handled exactly).12 S7 Newton basin of attraction in the exponential model: convergence iff 0 < x < 2 theta_hat.13 S8 Remainder-amplification constant A_r = sum |w_j| s_j and the sufficient exponent table alpha(r,q).14"""15import sympy as sp16from fractions import Fraction1718def vandermonde_weights(s):19 """Weights w solving sum w_j = 1, sum w_j s_j^a = 0 (a=1..r-1); exact rationals."""20 r = len(s)21 M = sp.Matrix([[sp.Integer(sj)**a for sj in s] for a in range(r)])22 rhs = sp.Matrix([1] + [0]*(r-1))23 return list(M.LUsolve(rhs))2425# ---------- S1: weights ----------26checks = {}27for s in [(1,2), (1,2,4), (1,2,4,8), (1,2,3), (1,3,9), (1,2,4,8,16)]:28 w = vandermonde_weights(s)29 for a in range(len(s)):30 tot = sum(wj * sp.Integer(sj)**a for wj, sj in zip(w, s))31 assert tot == (1 if a == 0 else 0), (s, a, tot)32 checks[s] = w33assert checks[(1,2)] == [2, -1]34assert checks[(1,2,4)] == [sp.Rational(8,3), -2, sp.Rational(1,3)]35print("S1 OK weights:", {k: [str(x) for x in v] for k, v in checks.items()})3637# ---------- S2: exponential model moments ----------38n, th = sp.symbols('n theta0', positive=True)39# S_n ~ Gamma(n, rate theta0): E[S^-k] = theta0^k * Gamma(n-k)/Gamma(n)40def Einv(k): # E[(theta0*Xbar_n)^(-k)] = E[(S_n/n)^-k] * theta0^-k ... compute E[Y^-k] with Y = theta0*Xbar_n, nY ~ Gamma(n,1)41 # E[(nY)^-k] = Gamma(n-k)/Gamma(n) => E[Y^-k] = n^k Gamma(n-k)/Gamma(n)42 return sp.simplify(n**k * sp.gamma(n-k)/sp.gamma(n))43EY1 = sp.simplify(sp.gammasimp(Einv(1)))44EY2 = sp.simplify(sp.gammasimp(Einv(2)))45assert sp.simplify(EY1 - n/(n-1)) == 046assert sp.simplify(EY2 - n**2/((n-1)*(n-2))) == 047# E[theta_hat_n] = theta0 * n/(n-1); Var = theta0^2 n^2/((n-1)^2 (n-2))48var = sp.simplify(th**2*(EY2 - EY1**2))49assert sp.simplify(var - th**2*n**2/((n-1)**2*(n-2))) == 050# remainder R_n = theta0 (Y^-1 - 2 + Y), E Y = 1, E Y^2 = 1 + 1/n51ER = sp.simplify(th*(EY1 - 2 + 1))52assert sp.simplify(ER - th/(n-1)) == 053ER2 = sp.simplify(th**2*(EY2 - 4*EY1 + 6 - 4*1 + (1 + 1/n)))54assert sp.simplify(ER2 - th**2*(3*n+2)/(n*(n-1)*(n-2))) == 055# series: E R_n = theta0 (1/n + 1/n^2 + ...) => b_a = theta0 for all a56ser = sp.series(th/(n-1), n, sp.oo, 4).removeO()57assert sp.simplify(ser - th*(1/n + 1/n**2 + 1/n**3)) == 058# uniform constant: sup_{n>=3} n^2 E R_n^2 / theta0^259f = sp.lambdify(n, n**2*(3*n+2)/(n*(n-1)*(n-2)))60Csup = max(f(k) for k in range(3, 2000))61assert abs(Csup - 16.5) < 1e-1262print("S2 OK E[theta_hat_n]=theta0 n/(n-1), Var=theta0^2 n^2/((n-1)^2(n-2)), E R_n=theta0/(n-1), E R_n^2=theta0^2(3n+2)/(n(n-1)(n-2)), sup n^2 E R_n^2/theta0^2 = 16.5 (n=3)")6364# ---------- S3: exact bias of J_r and leading coefficient ----------65m = sp.symbols('m', positive=True)66def exact_bias(s):67 w = vandermonde_weights(s)68 # sub-block size m/s_j, E[theta_hat] = theta0 (m/s_j)/((m/s_j)-1) = theta0 m/(m - s_j)69 return sp.simplify(sum(wj*th*m/(m-sj) for wj, sj in zip(w, s)) - th)70b2 = exact_bias((1,2)); b3 = exact_bias((1,2,4))71assert sp.simplify(b2 + 2*th/((m-1)*(m-2))) == 072assert sp.simplify(b3 - 8*th/((m-1)*(m-2)*(m-4))) == 073betas = {}74for s in [(1,2),(1,2,4),(1,2,4,8),(1,2,3),(1,3,9)]:75 w = vandermonde_weights(s); r = len(s)76 beta = sum(wj*sp.Integer(sj)**r for wj, sj in zip(w, s))77 lead = sp.limit(exact_bias(s)*m**r, m, sp.oo)78 assert sp.simplify(lead - th*beta) == 0, (s, lead, beta)79 betas[s] = beta80assert betas[(1,2)] == -2 and betas[(1,2,4)] == 8 and betas[(1,2,4,8)] == -6481print("S3 OK exact biases: J2 = -2 theta0/((m-1)(m-2)), J3 = 8 theta0/((m-1)(m-2)(m-4)); beta_r:", {k: str(v) for k, v in betas.items()})8283# ---------- S4: Newton identity ----------84x, xbar = sp.symbols('x xbar', positive=True)85L = sp.log(x) - x*xbar # average log-likelihood of Exp(rate x)86g = sp.diff(L, x); H = sp.diff(L, x, 2)87Nx = sp.simplify(x - g/H)88assert sp.simplify(Nx - (2*x - x**2*xbar)) == 089that = 1/xbar90assert sp.simplify((Nx - that) + (x - that)**2/that) == 091print("S4 OK N(x) = 2x - x^2 Xbar, N(x) - theta_hat = -(x - theta_hat)^2/theta_hat")9293# ---------- S5: q-fold identity and boundary limit ----------94e0, thh, c, N = sp.symbols('e0 thetahat c N', positive=True)95def newton_err(e, q):96 for _ in range(q):97 e = -e**2/thh98 return sp.simplify(e)99for q in range(1, 4):100 assert sp.simplify(newton_err(e0, q) + e0**(2**q)/thh**(2**q-1)) == 0101# boundary: K = c N^alpha, alpha = 1 - 1/(2^(q+1) r) => m = N/K = N^{1/(2^(q+1) r)}/c ; e0 ~ beta_r theta0 m^{-r}102for s in [(1,2),(1,2,4),(1,2,4,8)]:103 r = len(s); beta = betas[s]104 for q in range(0, 3):105 alpha = 1 - sp.Rational(1, 2**(q+1)*r)106 mm = N**(1-alpha)/c107 e_pilot = beta*th*mm**(-r)108 if q == 0:109 lim = sp.limit(sp.sqrt(N)*e_pilot * N**(sp.Rational(1,2)) / N**(sp.Rational(1,2)), N, sp.oo) # sqrt(N) e0 -> infinity: use N^{1/2} scaling only for q>=1110 err = newton_err(e_pilot, q).subs(thh, th)111 scaled = sp.simplify(sp.sqrt(N)*err)112 lim = sp.limit(scaled, N, sp.oo)113 expected = -th*(beta*c**r)**(2**q) if q >= 1 else sp.nan114 if q >= 1:115 assert sp.simplify(lim - expected) == 0, (s, q, lim, expected)116 else:117 # q = 0: the scaled bias N^{1/2} m^{-r} = c^r N^{1/2 - 1/2} -> beta theta0 c^r (alpha(r,0) = 1 - 1/(2r) makes sqrt(N) m^{-r} = c^r)118 assert sp.simplify(lim - beta*th*c**r) == 0, (s, q, lim)119assert sp.simplify(-th*((-2)*c**2)**2 + 4*th*c**4) == 0120print("S5 OK N^q(x) - theta_hat = -(x-theta_hat)^(2^q)/theta_hat^(2^q-1); boundary limit sqrt(N)(theta^(q) - theta_hat) -> -theta0 (beta_r c^r)^(2^q) for q>=1 (r=2,q=1: -4 theta0 c^4) and beta_r theta0 c^r for q=0")121122# ---------- S6: exact variance of J2 in the exponential model ----------123# S_A, S_B ~ Gamma(h, 1), h = m/2 (theta0 = 1 w.l.o.g.); S = S_A+S_B ~ Gamma(2h,1), W = S_A/S ~ Beta(h,h) indep. of S.124# theta_full = m/S = 2h/S ; theta_A = h/S_A = h/(S W) ; theta_B = h/(S (1-W)).125# J2 = 2*theta_full - (theta_A+theta_B)/2 = 4h/S - (h/2)(1/W + 1/(1-W))/S126h = sp.symbols('h', positive=True)127def EgammaInv(a, k): # E[S^-k], S ~ Gamma(a,1)128 return sp.gamma(a-k)/sp.gamma(a)129def EbetaInv(a, b, k, l): # E[W^-k (1-W)^-l], W ~ Beta(a,b)130 return sp.gamma(a-k)*sp.gamma(b-l)*sp.gamma(a+b)/(sp.gamma(a)*sp.gamma(b)*sp.gamma(a+b-k-l))131# J2 = (1/S) * [4h - (h/2) (1/W + 1/(1-W))] =: (1/S) * G(W)132# E J2 = E[S^-1] * E[G]; E J2^2 = E[S^-2] * E[G^2]133EG = 4*h - (h/2)*(EbetaInv(h,h,1,0) + EbetaInv(h,h,0,1))134EG2 = 16*h**2 - 2*4*h*(h/2)*(EbetaInv(h,h,1,0)+EbetaInv(h,h,0,1)) + (h**2/4)*(EbetaInv(h,h,2,0) + EbetaInv(h,h,0,2) + 2*EbetaInv(h,h,1,1))135EJ2 = sp.gammasimp(EgammaInv(2*h,1)*EG)136EJ22 = sp.gammasimp(EgammaInv(2*h,2)*EG2)137EJ2_m = sp.simplify(EJ2.subs(h, m/2))138VarJ2_m = sp.simplify((EJ22 - EJ2**2).subs(h, m/2))139assert sp.simplify(EJ2_m - (1 - 2/((m-1)*(m-2)))) == 0 # matches S3 (theta0 = 1)140# asymptotic expansion of the variance: 1/m + v2/m^2 + ...141serV = sp.series(VarJ2_m, m, sp.oo, 4).removeO()142v1 = sp.limit(VarJ2_m*m, m, sp.oo); v2 = sp.limit((VarJ2_m - 1/m)*m**2, m, sp.oo)143assert v1 == 1144# compare with the full-block MLE variance m^2/((m-1)^2 (m-2)) = 1/m + 3/m^2 + ...145v2_full = sp.limit((m**2/((m-1)**2*(m-2)) - 1/m)*m**2, m, sp.oo)146print("S6 OK Var(J2) exact =", sp.factor(VarJ2_m), " ~ 1/m + (%s)/m^2 ; Var(theta_hat_m) ~ 1/m + (%s)/m^2" % (v2, v2_full))147VARJ2_EXACT = VarJ2_m148149# ---------- S7: Newton basin ----------150# e_{j+1} = -e_j^2/theta_hat ; |e_{j+1}|/theta_hat = (|e_j|/theta_hat)^2 -> 0 iff |e_0| < theta_hat iff 0 < x < 2 theta_hat151u = sp.symbols('u', positive=True)152assert sp.simplify(sp.Abs(-(u*thh)**2/thh)/thh - u**2) == 0153print("S7 OK Newton basin: |x - theta_hat| < theta_hat <=> 0 < x < 2 theta_hat ; x <= 0 or x >= 2 theta_hat diverges (x = 2 theta_hat maps to 0 then stays at 0)")154155# ---------- S8: amplification constants and exponent table ----------156def A_const(s):157 w = vandermonde_weights(s); return sum(abs(wj)*sj for wj, sj in zip(w, s))158def W_sum_abs(s):159 w = vandermonde_weights(s); return sum(abs(wj) for wj in w)160tab = []161for s in [(1,), (1,2), (1,2,4), (1,2,4,8), (1,2,3), (1,3,9)]:162 tab.append((s, A_const(s), W_sum_abs(s)))163assert A_const((1,2)) == 4 and A_const((1,2,4)) == 8164alpha = {(r,q): sp.Rational(1) - sp.Rational(1, 2**(q+1)*r) for r in (1,2,3,4) for q in (0,1,2,3)}165assert alpha[(1,0)] == sp.Rational(1,2) and alpha[(1,1)] == sp.Rational(3,4) and alpha[(2,1)] == sp.Rational(7,8) and alpha[(3,1)] == sp.Rational(11,12) and alpha[(1,2)] == sp.Rational(7,8) and alpha[(2,0)] == sp.Rational(3,4) and alpha[(3,0)] == sp.Rational(5,6)166print("S8 OK A_r = sum|w_j| s_j:", [(s, str(a), str(b)) for s, a, b in tab])167print(" alpha(r,q) = 1 - 1/(2^(q+1) r):", {k: str(v) for k, v in alpha.items()})168print("ALL SYMBOLIC CHECKS PASSED")169