Statistical Validation Suite: AIC/BIC, Leave-One-Out CV, Bootstrap, Permutation Test, and Moran's I Spatial Autocorrelation

#!/usr/bin/env python3
"""
Statistical Validation Suite — Mildner (2026) Germania Magna v8
Kinematic Block-Deformation Model
================================================================
Implements:
  1. Formal AIC, AICc, BIC (Akaike 1974; Burnham & Anderson 2002)
  2. Wilcoxon signed-rank test on blind-test residuals
  3. Leave-One-Out Cross-Validation (LOO-CV) for EC cluster
  4. Bootstrap confidence intervals for delta-lambda_EC
  5. Permutation test for cluster spatial specificity
  6. Monte Carlo robustness under identification uncertainty
Data: Mildner (2026a,b); Tables 1 & 4
      EarthArXiv DOI: 10.31223/X5KB51 / 10.31223/X5313T
      Zenodo: 10.5281/zenodo.10968193
Requirements: numpy, scipy, matplotlib
Usage:        python statistical_validation_suite.py
================================================================
"""
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import warnings
warnings.filterwarnings('ignore')
np.random.seed(2026)   # reproducible
# ── DATA ──────────────────────────────────────────────────────
# 22 non-calibration evaluation points (v8)
# r_B for training points: estimated from known corrections;
# r_B for 7 blind-test points: exact (Table 4).
# S-C: v8 Bernsdorf identification (same Δλ = -70.0 km as v7.3)
# S6/S7: R_B_BERNSDORF override values (v8 primary)
# Columns: name, Δλ_A (km), Δφ_A (km), r_A (km), r_B (km),
#          in_EC_cluster, is_blind_test_point
POINTS = [
    # Training points (15)
    ('S3_Bud',   -87.1,  26.9,  91.2,  34.0, True,  False),
    ('S4_Cal',   -38.2,  13.7,  40.6,  38.5, False, False),
    ('SA_Ars',   -51.5,  17.0,  54.2,  52.0, False, False),
    ('S5_Lim',   -78.2,   9.2,  78.7,  28.0, True,  False),
    ('SL_Leu',   -87.4,  23.7,  90.5,  32.0, True,  False),
    ('SC_Car',   -70.0,   3.6,  70.1,  22.0, True,  False),   # v8 Bernsdorf
    ('S8_Tre',   +26.6,  -7.9,  27.8,  27.8, False, False),
    ('S9_Lir',    -5.9, -24.6,  25.3,  25.3, False, False),
    ('G1_ANW',   -46.5,  17.6,  49.7,  45.0, False, False),
    ('G4_MeE',   -41.2,  39.1,  55.2,  42.0, False, False),
    ('G5_SuW',   -19.8, -13.7,  23.9,  11.0, False, False),
    ('G7_Sar',   -93.9,  -2.2,  94.0,  90.0, False, False),
    ('G8_Abn',   +11.2,  91.2,  91.9,  11.2, False, False),
    ('F4_Sue',   -63.4,  16.4,  65.5,  60.0, False, False),
    ('F5_Via',   -32.3,   7.1,  33.1,  31.0, False, False),
    # Blind-test points (7) — v8 primary values
    ('S6_Lug',  -109.5,  24.5, 112.2,  31.8, True,  True),   # v8 Bernsdorf r_B
    ('S7_Str',  -101.0,  11.8, 101.7,  35.2, True,  True),   # v8 Bernsdorf r_B
    ('G6_SuE',   +11.2,  71.6,  72.5,   7.5, False, True),
    ('G3_MeW',   -48.6,  -1.3,  48.6,  62.5, False, True),
    ('G2_ASE',   -31.3,  18.7,  36.3,  36.8, False, True),
    ('F3_Cha',   -74.5,  21.0,  77.4,  77.4, False, True),
    ('F2_Vis',  -124.0,  64.0, 142.0, 127.2, False, True),
]
names  = [p[0] for p in POINTS]
dL_A   = np.array([p[1] for p in POINTS])
r_A    = np.array([p[3] for p in POINTS])
r_B    = np.array([p[4] for p in POINTS])
EC     = np.array([p[5] for p in POINTS], dtype=bool)
test   = np.array([p[6] for p in POINTS], dtype=bool)
n     = len(POINTS)
k_A   = 6
k_B   = 9
EC_dL      = dL_A[EC]
EC_names   = np.array(names)[EC]
r_A_test   = r_A[test]
r_B_test   = r_B[test]
test_names = np.array(names)[test]

# ── HELPER ────────────────────────────────────────────────────
def info_criteria(residuals, k, n_obs):
    RSS   = np.sum(residuals**2)
    nlogL = n_obs * np.log(RSS / n_obs)
    AIC   = nlogL + 2 * k
    denom = n_obs - k - 1
    AICc  = AIC + (2 * k * (k + 1) / denom) if denom > 0 else np.inf
    BIC   = nlogL + np.log(n_obs) * k
    return dict(AIC=AIC, AICc=AICc, BIC=BIC, RSS=RSS, RMSE=np.sqrt(RSS / n_obs))
def evidence_label(d):
    if d < 2:  return "inconclusive"
    if d < 6:  return "moderate support for B"
    if d < 10: return "strong support for B"
    return            "VERY STRONG support for B"

# ── PART 1: AIC / AICc / BIC ──────────────────────────────────
mA = info_criteria(r_A, k_A, n)
mB = info_criteria(r_B, k_B, n)
dAIC  = mA['AIC']  - mB['AIC']
dAICc = mA['AICc'] - mB['AICc']
dBIC  = mA['BIC']  - mB['BIC']
mB_k8  = info_criteria(r_B, 8, n)
dAICc8 = mA['AICc'] - mB_k8['AICc']
w_raw_A = np.exp(-0.5 * (mA['AIC'] - min(mA['AIC'], mB['AIC'])))
w_raw_B = np.exp(-0.5 * (mB['AIC'] - min(mA['AIC'], mB['AIC'])))
wA = w_raw_A / (w_raw_A + w_raw_B)
wB = w_raw_B / (w_raw_A + w_raw_B)

# ── PART 2: WILCOXON ──────────────────────────────────────────
diffs = r_A_test - r_B_test
stat_w, p_w = stats.wilcoxon(diffs, alternative='greater')
r_B_corr = r_B_test.copy()
g3 = list(test_names).index('G3_MeW')
r_B_corr[g3] = 48.6
diffs_c = r_A_test - r_B_corr
stat_wc, p_wc = stats.wilcoxon(diffs_c, alternative='greater')

# ── PART 3: LOO-CV (EC CLUSTER) ───────────────────────────────
loo_pred = np.array([EC_dL[np.arange(len(EC_dL)) != i].mean()
                     for i in range(len(EC_dL))])
loo_err  = EC_dL - loo_pred
loo_rmse = np.sqrt(np.mean(loo_err**2))
null_rmse_loo = np.sqrt(np.mean(EC_dL**2))

# ── PART 4: BOOTSTRAP CI ──────────────────────────────────────
N_BOOT = 200_000
boot_means = np.array([
    np.random.choice(EC_dL, len(EC_dL), replace=True).mean()
    for _ in range(N_BOOT)
])
ci95   = np.percentile(boot_means, [2.5, 97.5])
ci99   = np.percentile(boot_means, [0.5, 99.5])
p_zero = np.mean(boot_means >= 0)

# ── PART 5: PERMUTATION TEST ───────────────────────────────────
N_PERM   = 500_000
n_EC     = EC.sum()
obs_mean = EC_dL.mean()
perm_means = np.array([
    np.random.choice(dL_A, n_EC, replace=False).mean()
    for _ in range(N_PERM)
])
p_perm = np.mean(perm_means <= obs_mean)
z_perm = (obs_mean - perm_means.mean()) / perm_means.std()

# ── PART 6: MONTE CARLO ROBUSTNESS ────────────────────────────
sigma_id  = 15.0
N_MC      = 100_000
t_crit001 = stats.t.ppf(0.0005, df=n_EC - 1)
t_crit005 = stats.t.ppf(0.025,  df=n_EC - 1)
mc_t = np.array([
    (lambda x: x.mean() / (x.std(ddof=1) / np.sqrt(n_EC)))(
        EC_dL + np.random.normal(0, sigma_id, n_EC))
    for _ in range(N_MC)
])
f001 = np.mean(mc_t < t_crit001)
f005 = np.mean(mc_t < t_crit005)

# ── FIGURE (six panels) ───────────────────────────────────────
CA, CB = '#E53935', '#1E88E5'
fig = plt.figure(figsize=(18, 12))
gs  = gridspec.GridSpec(2, 3, figure=fig, hspace=0.45, wspace=0.35)
# [Panel code as in v7.3; observed t updated to -15.0 for v8]
fig.suptitle(
    "Statistical Validation Suite — Mildner (2026) Kinematic Germania Magna Model v8\n"
    "AIC/BIC  ·  Bootstrap  ·  Permutation Test  ·  LOO-CV  ·  Monte Carlo Robustness",
    fontsize=12, fontweight='bold')
plt.savefig('mildner_stat_validation.png', dpi=150,
            bbox_inches='tight', facecolor='white')
print("Figure saved: mildner_stat_validation.png")

#!/usr/bin/env python3
"""
Moran's I Spatial Autocorrelation — Mildner (2026) Germania Magna v8
Falsification Test T39
================================================================
Tests whether cartometric residuals of Ptolemy's Germania Magna
(Model B, v8) are spatially autocorrelated, as predicted by the
block-displacement hypothesis.
  H0 (uniform-error null):  residuals are i.i.d.; E[I] = -1/(n-1)
  H1 (block-displacement):  residuals are positively autocorrelated
Weight schemes tested:
  (1) 1/d²  (raw)
  (2) 1/d   (raw)
  (3) 1/d²  (row-standardised)
  (4) kNN k=4 (row-standardised)
Variable x_i:
  (a) residual magnitude r_i
  (b) longitudinal residual Δλ_i
Produces: morans_i_diagnostics.png  (four panels)
Reference: Appendix D, Mildner (2026b), EarthArXiv v8
           Zenodo: 10.5281/zenodo.10968193
Requirements: numpy, scipy, matplotlib
Usage:        python morans_i_test.py
================================================================
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import warnings
warnings.filterwarnings('ignore')
np.random.seed(2026)
# ── DATA: 22-point v8 dataset (Table D_dataset, Appendix D) ──
# Columns: ID, lon (°E), lat (°N), Δλ (km), r_magnitude (km), EC
POINTS = [
    ('F1',   14.33, 51.31,  -61.5,  80.5,  False),
    ('F2',   13.82, 51.20, -124.0, 142.0,  False),
    ('F3',   12.08, 52.83,  -74.5,  77.4,  False),
    ('F4',   12.81, 52.92,  -63.4,  65.5,  False),
    ('F5',   13.68, 52.85,  -32.3,  33.1,  False),
    ('S3',   13.56, 51.62,  -87.1,  91.2,  True),
    ('S4',   13.99, 51.74,  -38.2,  40.6,  False),
    ('SA',   14.00, 51.52,  -51.5,  54.2,  False),
    ('S5',   13.50, 52.05,  -78.2,  78.7,  True),
    ('S6',   13.24, 51.58, -109.5, 112.2,  True),
    ('S7',   13.23, 51.69, -101.0, 101.7,  True),
    ('SL',   13.71, 51.63,  -87.4,  90.5,  True),
    ('SC',   14.05, 51.38,  -70.0,  70.0,  True),   # v8 Bernsdorf
    ('S8',    8.80, 53.08,  +26.6,  27.8,  False),
    ('S9',    9.95, 53.06,   -5.9,  25.3,  False),
    ('G1',   11.70, 52.13,  -46.5,  49.7,  False),
    ('G2',   13.95, 51.77,  -31.3,  36.3,  False),
    ('G3',    9.83, 51.99,  -48.6,  48.6,  False),
    ('G4',   11.55, 51.53,  -41.2,  55.2,  False),
    ('G5',    9.49, 51.31,  -19.8,  23.9,  False),
    ('G6',   11.00, 50.70,  +11.2,  72.5,  False),
    ('G7',   14.50, 51.00,  -93.9,  94.0,  False),
]
ids     = [p[0] for p in POINTS]
lons    = np.array([p[1] for p in POINTS])
lats    = np.array([p[2] for p in POINTS])
dl      = np.array([p[3] for p in POINTS])    # Δλ (km)
rmag    = np.array([p[4] for p in POINTS])    # residual magnitude (km)
ec_mask = np.array([p[5] for p in POINTS])
n       = len(POINTS)

# ── GREAT-CIRCLE DISTANCE MATRIX (km) ────────────────────────
def great_circle_km(lons_deg, lats_deg):
    R    = 6371.0
    lo   = np.radians(lons_deg)
    la   = np.radians(lats_deg)
    n_pt = len(lo)
    D    = np.zeros((n_pt, n_pt))
    for i in range(n_pt):
        for j in range(i + 1, n_pt):
            dlat = la[j] - la[i]
            dlon = lo[j] - lo[i]
            a = (np.sin(dlat / 2) ** 2
                 + np.cos(la[i]) * np.cos(la[j]) * np.sin(dlon / 2) ** 2)
            D[i, j] = D[j, i] = 2 * R * np.arcsin(np.sqrt(a))
    return D
D = great_circle_km(lons, lats)

# ── WEIGHT MATRICES ──────────────────────────────────────────
def build_weights(D_mat, scheme, k=4):
    n_pt = len(D_mat)
    W    = np.zeros((n_pt, n_pt))
    with np.errstate(divide='ignore', invalid='ignore'):
        if scheme == 'inv_d2':
            W = np.where(D_mat > 0, 1.0 / D_mat ** 2, 0.0)
        elif scheme == 'inv_d':
            W = np.where(D_mat > 0, 1.0 / D_mat, 0.0)
        elif scheme == 'inv_d2_rs':
            raw      = np.where(D_mat > 0, 1.0 / D_mat ** 2, 0.0)
            row_sums = raw.sum(axis=1, keepdims=True)
            W        = np.where(row_sums > 0, raw / row_sums, 0.0)
        elif scheme == 'knn_rs':
            raw = np.where(D_mat > 0, 1.0 / D_mat ** 2, 0.0)
            for i in range(n_pt):
                row  = raw[i].copy()
                thr  = np.sort(row)[::-1][k] if n_pt > k else 0
                row[row < thr] = 0.0
                W[i] = row
            row_sums = W.sum(axis=1, keepdims=True)
            W        = np.where(row_sums > 0, W / row_sums, 0.0)
    np.fill_diagonal(W, 0.0)
    return W

# ── MORAN'S I ────────────────────────────────────────────────
def morans_i(x, W):
    """Return (I, E_I) for variable x and weight matrix W."""
    n_pt  = len(x)
    xc    = x - x.mean()
    W_sum = W.sum()
    num   = float(np.einsum('ij,i,j->', W, xc, xc))   # Σ_ij w_ij x_i x_j
    den   = float(xc @ xc)
    I     = (n_pt / W_sum) * (num / den)
    E_I   = -1.0 / (n_pt - 1)
    return I, E_I

def analytic_z(x, W_mat, I_obs):
    """Analytic z-score under randomisation assumption."""
    n_pt  = len(x)
    E_I   = -1.0 / (n_pt - 1)
    S1    = 0.5 * float(np.sum((W_mat + W_mat.T) ** 2))
    S2    = float(np.sum((W_mat.sum(axis=1) + W_mat.sum(axis=0)) ** 2))
    W_sum = W_mat.sum()
    n2    = n_pt ** 2
    m2    = float(np.sum((x - x.mean()) ** 2)) / n_pt
    m4    = float(np.sum((x - x.mean()) ** 4)) / n_pt
    b2    = m4 / m2 ** 2
    v_num = n_pt * ((n2 - 3 * n_pt + 3) * S1 - n_pt * S2 + 3 * W_sum ** 2)
    v_den = ((n_pt - 1)
             * ((n2 + 3 * n_pt - 6) * S1
                - (n2 - n_pt + 2) * S2
                + 6 * W_sum ** 2))
    var_I = (n_pt * v_num - b2 * v_den) / ((n_pt - 1) ** 2 * (n_pt + 1) * W_sum ** 2)
    if var_I <= 0:
        return 0.0
    return (I_obs - E_I) / np.sqrt(var_I)

def permutation_p(x, W_mat, I_obs, n_perm=9999):
    """One-sided permutation p-value P(I_rand >= I_obs)."""
    count = 0
    xp    = x.copy()
    for _ in range(n_perm):
        np.random.shuffle(xp)
        I_rand, _ = morans_i(xp, W_mat)
        if I_rand >= I_obs:
            count += 1
    return (count + 1) / (n_perm + 1)

# ── RUN TESTS ────────────────────────────────────────────────
SCHEMES = [
    ('inv_d2',    '1/d² (raw)'),
    ('inv_d',     '1/d (raw)'),
    ('inv_d2_rs', '1/d² row-stand.'),
    ('knn_rs',    'kNN k=4 row-stand.'),
]
E_I_null = -1.0 / (n - 1)
print("=" * 72)
print("MORAN'S I — Falsification Test T39")
print("Mildner (2026) Germania Magna v8  |  n=22  seed=2026")
print("=" * 72)
print(f"\nE[I] under H0 = {E_I_null:.3f}\n")
hdr = (f"{'Scheme':<22}  {'I(r)':>6}  {'z(r)':>6}  {'p(r)':>7}  "
       f"{'I(Δλ)':>6}  {'z(Δλ)':>6}  {'p(Δλ)':>7}")
print(hdr)
print("-" * 72)
results = {}
for scheme_key, scheme_label in SCHEMES:
    W_mat   = build_weights(D, scheme_key)
    I_r,  _ = morans_i(rmag, W_mat)
    I_dl, _ = morans_i(dl,   W_mat)
    z_r     = analytic_z(rmag, W_mat, I_r)
    z_dl    = analytic_z(dl,   W_mat, I_dl)
    p_r     = permutation_p(rmag, W_mat, I_r)
    p_dl    = permutation_p(dl,   W_mat, I_dl)
    sig_r   = '*' if p_r  < 0.05 else ' '
    sig_dl  = '*' if p_dl < 0.05 else ' '
    print(f"{scheme_label:<22}  {I_r:>+6.3f}  {z_r:>+6.2f}  "
          f"{p_r:>6.3f}{sig_r}  {I_dl:>+6.3f}  {z_dl:>+6.2f}  "
          f"{p_dl:>6.3f}{sig_dl}")
    results[scheme_key] = (I_r, z_r, p_r, I_dl, z_dl, p_dl)
print("-" * 72)
print("* p < 0.05 (one-sided permutation, 9999 shuffles, seed=2026)")
n_sig_r  = sum(1 for v in results.values() if v[2] < 0.05)
n_sig_dl = sum(1 for v in results.values() if v[5] < 0.05)
print(f"\nH0 rejected at p<0.05: {n_sig_r}/4 schemes (r)  "
      f"and {n_sig_dl}/4 (Δλ).")
print("\nStatus: PRELIMINARY — peer-reviewed coordinates and")
print("pre-registered weight scheme required for T39 confirmatory result.")

# ── FIGURE: four diagnostic panels ───────────────────────────
fig = plt.figure(figsize=(14, 10))
gs  = gridspec.GridSpec(2, 2, figure=fig, hspace=0.42, wspace=0.38)
# (a) Spatial distribution of residual magnitudes
ax0 = fig.add_subplot(gs[0, 0])
sc  = ax0.scatter(lons[~ec_mask], lats[~ec_mask],
                  c=rmag[~ec_mask], s=60, cmap='YlOrRd',
                  vmin=20, vmax=145, edgecolors='k', lw=0.5)
ax0.scatter(lons[ec_mask], lats[ec_mask],
            c=rmag[ec_mask], s=100, cmap='YlOrRd',
            vmin=20, vmax=145, edgecolors='blue', lw=1.5,
            marker='D', label='Elster Cluster (EC)')
plt.colorbar(sc, ax=ax0, label='Residual magnitude $r_i$ (km)')
ax0.set_xlabel('Longitude (°E)')
ax0.set_ylabel('Latitude (°N)')
ax0.set_title('(a) Spatial distribution of $r_i$')
ax0.legend(fontsize=8)
# (b,c) Moran scatterplots — row-standardised 1/d²
W_rs = build_weights(D, 'inv_d2_rs')
for ax_i, (x_var, var_label, panel) in enumerate([
    (rmag, '$r_i$ (km)',           '(b)'),
    (dl,   r'$\Delta {\lambda }_i$ (km)', '(c)'),
]):
    ax = fig.add_subplot(gs[0, 1] if ax_i == 0 else gs[1, 0])
    Wx = W_rs @ x_var
    xc = x_var - x_var.mean()
    m  = np.polyfit(xc, Wx - x_var.mean(), 1)
    ax.scatter(xc, Wx - x_var.mean(), c='steelblue', s=40, alpha=0.75)
    ax.plot(np.sort(xc), np.polyval(m, np.sort(xc)),
            'r-', lw=1.5, label=f'slope = I = {m[0]:.3f}')
    ax.axhline(0, lw=0.6, color='grey')
    ax.axvline(0, lw=0.6, color='grey')
    ax.set_xlabel(f'Deviation of {var_label}')
    ax.set_ylabel(f'Spatial lag of {var_label}')
    ax.set_title(f'{panel} Moran scatterplot — row-std 1/d²\n{var_label}')
    ax.legend(fontsize=8)
# (d) Permutation reference distributions for all four schemes
ax3   = fig.add_subplot(gs[1, 1])
cols  = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']
xp_r  = rmag.copy()
for (sk, sl), col in zip(SCHEMES, cols):
    W_t       = build_weights(D, sk)
    I_obs_t,_ = morans_i(rmag, W_t)
    dist = []
    for _ in range(999):
        np.random.shuffle(xp_r)
        Ip, _ = morans_i(xp_r, W_t)
        dist.append(Ip)
    ax3.hist(dist, bins=30, alpha=0.35, color=col, label=sl)
    ax3.axvline(I_obs_t, color=col, lw=2, linestyle='--')
ax3.set_xlabel("Moran's $I$")
ax3.set_ylabel('Frequency (permutation, 999 shuffles)')
ax3.set_title("(d) Permutation distributions\n"
              "dashed lines = observed $I$ per scheme")
ax3.legend(fontsize=7)
fig.suptitle(
    f"Moran's I Spatial Autocorrelation — Germania Magna v8 (T39)\n"
    f"n={n},  E[I]={E_I_null:.3f},  seed=2026  (preliminary)",
    fontsize=11)
plt.savefig('morans_i_diagnostics.png', dpi=150, bbox_inches='tight')
print("\nFigure saved: morans_i_diagnostics.png")

=================================================================
PART 1  INFORMATION CRITERIA
=================================================================
Model A — Affine (k=6)
  RMSE = 73.98 km   RSS = 120393 km²
  AIC  = 201.36   AICc = 206.96   BIC = 207.91
Model B — Kinematic (k=9)
  RMSE = 50.07 km   RSS = 55145 km²
  AIC  = 190.19   AICc = 205.19   BIC = 200.01
  ΔAIC  = +11.18  →  VERY STRONG support for B
  ΔAICc = +1.78   →  inconclusive
  ΔAICc (c as prior, k_eff=8) = +7.70  →  strong support for B
  ΔBIC  = +7.90   →  strong support for B
  Akaike weights: w(A)=0.0037  w(B)=0.9963
  Model B is 267× more probable (AIC)
  ⚠  AICc note: The standard AICc penalises all k=9 parameters
     against n=22. However:
     • c = 15.2 km/°P is a geological prior (Taunus stability),
       not estimated from the 22-point dataset.
     • δλ_EC was estimated from 4 held-in training points only.
     • θ was estimated from G5 (1 training point) alone.
     Using k_eff=8 (c as prior) gives ΔAICc = +7.7 (strong
     support). The blind-test and permutation results below are
     independent of this parameter-counting debate.
=================================================================
PART 2  WILCOXON SIGNED-RANK TEST — BLIND TEST (n=7, v8)
=================================================================
  Point          r_A     r_B(v8)     Δ
  S6_Lug       112.2      31.8    +80.4  ▲
  S7_Str       101.7      35.2    +66.5  ▲
  G6_SuE        72.5       7.5    +65.0  ▲
  G3_MeW        48.6      62.5    -13.9  ▼
  G2_ASE        36.3      36.8     -0.5  ▼
  F3_Cha        77.4      77.4     +0.0  =
  F2_Vis       142.0     127.2    +14.8  ▲
  RMSE Model A: 91.0 km  →  Model B (v8): 65.2 km  (-28.4%)
  [v7.3 reference: RMSE_B = 64.8 km  (-28.8%)]
  Wilcoxon W = 18, p = 0.0781
  With c_stable=0 for G3 (physically justified):
  Wilcoxon W = 14, p = 0.0625  (marginal)
=================================================================
PART 3  LEAVE-ONE-OUT CV — EC CLUSTER TRANSLATION (v8 Bernsdorf)
=================================================================
  Point        Observed   LOO-pred    Error
  S3_Bud          -87.1      -89.2     +2.1
  S5_Lim          -78.2      -91.0    +12.8
  SL_Leu          -87.4      -89.2     +1.8
  SC_Car          -70.0      -92.6    +22.6   ← v8 Bernsdorf
  S6_Lug         -109.5      -84.7    -24.8
  S7_Str         -101.0      -86.4    -14.6
  LOO-CV RMSE:  15.9 km
  Null RMSE:    89.8 km  (predict Δλ = 0)
  Improvement over null: 82%
  Mean bias:    0.00 km  (translation estimate: unbiased)
=================================================================
PART 4  BOOTSTRAP CI — δλ_EC  (n_boot = 200 000)
=================================================================
  Observed mean  = -88.87 km   SD = 14.48 km
  Bootstrap SE   = 5.41 km
  95% CI:  [-99.3, -78.5] km
  99% CI:  [-102.9, -75.6] km
  P(δλ_EC ≥ 0):  0.00e+00  (zero rigorously excluded)
=================================================================
PART 5  PERMUTATION TEST — EC CLUSTER SPECIFICITY
=================================================================
  H₀: The 6 EC points are a random sample of all 22 points
  Observed EC mean Δλ:  -88.9 km
  Population mean Δλ:   -52.5 km
  Permutation mean:     -52.5 km  SD = 14.3 km
  z = -2.55   p = 0.0030  (one-tailed)
=================================================================
PART 6  MONTE CARLO — ID UNCERTAINTY ROBUSTNESS
=================================================================
  ID uncertainty: ±15 km  (conservative upper bound)
  n_MC = 100,000
  EC cluster t (v8 observed): -15.0
  MC mean t = -11.99
  Runs achieving p<0.001: 97%
  Runs achieving p<0.05:  100%
=================================================================
SUMMARY TABLE
=================================================================
Test                          Result              Interpretation
─────────────────────────────────────────────────────────────────────
AIC  (n=22, k: 6 vs 9)       ΔAIC  = +11.2      Very strong (>10)
AICc (standard)               ΔAICc = +1.8       Inconclusive
AICc (c as prior, k=8)       ΔAICc = +7.7       Strong (>6)
BIC                            ΔBIC  = +7.9       Strong (>6)
Wilcoxon (7 test pts, v8)     p = 0.078          Uninformative (n=7)
Wilcoxon (G3 corrected)       p = 0.063          Uninformative (n=7)
EC LOO-CV RMSE                16 km vs 90 km     82% reduction
Bootstrap 95% CI              [-99, -78] km      Excludes zero
Permutation test               p = 0.0030         EC cluster non-random
MC ±15 km robustness          97% runs p<0.001   Robust
─────────────────────────────────────────────────────────────────────

========================================================================
MORAN'S I — Falsification Test T39
Mildner (2026) Germania Magna v8  |  n=22  seed=2026
========================================================================
E[I] under H0 = -0.048
Scheme                    I(r)    z(r)    p(r)    I(Δλ)   z(Δλ)   p(Δλ)
------------------------------------------------------------------------
1/d² (raw)              +0.586   +1.42   0.065   +0.361   +0.91   0.130
1/d (raw)               +0.184   +2.51   0.011*  +0.188   +2.55   0.013*
1/d² row-stand.         +0.352   +3.28   0.001*  +0.349   +3.24   0.001*
kNN k=4 row-stand.      +0.269   +2.69   0.012*  +0.321   +3.13   0.005*
------------------------------------------------------------------------
* p < 0.05 (one-sided permutation, 9999 shuffles, seed=2026)
H0 rejected at p<0.05: 3/4 schemes (r) and 3/4 (Δλ).
Status: PRELIMINARY — peer-reviewed coordinates and
pre-registered weight scheme required for T39 confirmatory result.
Figure saved: morans_i_diagnostics.png