Experiment 013 — Resampling Convergence Analysis

Studies how quickly bootstrap and RAWR confidence intervals converge as the number of replicates increases, to answer:

1. How many replicates are needed for a stable CI?
2. Does RAWR converge faster or slower than standard bootstrap?
3. How does data distribution affect convergence rate?
4. Is there a diminishing-returns threshold?

Method:

• Run bootstrap and RAWR at geometrically increasing replicate counts
• Track CI width convergence
• Measure coverage at each replicate count
• Compare convergence across Gaussian, log-normal, and heavy-tailed data

Reference: Lee & Liu (2024) IEEE BIBM — statistical resampling optimization Wang et al. (2021) Bioinformatics (ISMB) 37:i111-i119

Cross-spring: upgrades MC methodology for all experiments.

Domain: Statistics Reference: CLT convergence for bootstrap resampling

Data source: control/resampling_convergence/resampling_convergence.py + benchmark_*.json

This notebook is the publication-grade Python baseline for Experiment 013. The identical computations are validated in Rust (see validate_* binary) and delegated to barraCuda for GPU acceleration.

import json
import math
import sys
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt

# Wire path to groundSpring control/ for common utilities
CONTROL = Path('..') / '..' / 'control'
sys.path.insert(0, str(CONTROL))
from common import *  # noqa: F403 — validation harness

# Load benchmark data
benchmark_path = CONTROL / 'resampling_convergence' / 'benchmark_resampling_convergence.json'
with open(benchmark_path) as f:
    benchmark = json.load(f)

PASS_COLOR = '#2ecc71'
FAIL_COLOR = '#e74c3c'
INFO_COLOR = '#3498db'

print(f'Loaded benchmark: benchmark_resampling_convergence.json')
print(f'Provenance: {benchmark.get("_provenance", {})}')

Validation

Initialization

model = benchmark["model"]
exp = benchmark["expected_results"]
replicate_counts = model["replicate_counts"]
confidence = model["confidence"]
data_n = model["data_n"]

print("groundSpring Exp 013: Resampling Convergence (Lee & Liu 2024)")
print(f"  Replicate counts: {replicate_counts}")
print(f"  Data size: {data_n}, Confidence: {confidence}")
print("  Cross-spring: all springs (MC methodology optimization)")

Gaussian (μ=5.0, σ=2.0

gauss = model["gaussian"]
rng_data = np.random.default_rng(gauss["seed"])
data_gauss = rng_data.normal(gauss["mu"], gauss["sigma"], data_n)

boot_widths_g = []
rawr_widths_g = []
for n_boot in replicate_counts:
    rng_b = np.random.default_rng(gauss["seed"] + n_boot)
    rng_r = np.random.default_rng(gauss["seed"] + n_boot + 50000)
    bw = bootstrap_ci_width(data_gauss, n_boot, confidence, rng_b)
    rw = rawr_ci_width(data_gauss, n_boot, confidence, rng_r)
    boot_widths_g.append(bw)
    rawr_widths_g.append(rw)
    print(f"  n={n_boot:5d}: bootstrap={bw:.4f}  RAWR={rw:.4f}")

check_true(
    "Bootstrap width decreasing (Gaussian)",
    boot_widths_g[-1] <= boot_widths_g[0] * 1.1,
)
check_true(
    "RAWR width decreasing (Gaussian)",
    rawr_widths_g[-1] <= rawr_widths_g[0] * 1.1,
)

rel_change_boot = abs(boot_widths_g[-1] - boot_widths_g[-2]) / max(boot_widths_g[-2], 1e-10)
rel_change_rawr = abs(rawr_widths_g[-1] - rawr_widths_g[-2]) / max(rawr_widths_g[-2], 1e-10)
print(f"  Relative change 5k→10k: bootstrap={rel_change_boot:.4f} RAWR={rel_change_rawr:.4f}")

check_max(
    "Bootstrap converged (5k→10k < 15%)",
    rel_change_boot, exp["relative_width_change_5k_to_10k_max"],
)
check_max(
    "RAWR converged (5k→10k < 15%)",
    rel_change_rawr, exp["relative_width_change_5k_to_10k_max"],
)

Log-Normal (μ_ln=1.0, σ_ln=0.8

lognorm = model["lognormal"]
rng_data2 = np.random.default_rng(lognorm["seed"])
data_ln = rng_data2.lognormal(lognorm["mu_ln"], lognorm["sigma_ln"], data_n)

boot_widths_ln = []
for n_boot in replicate_counts:
    rng_b = np.random.default_rng(lognorm["seed"] + n_boot)
    bw = bootstrap_ci_width(data_ln, n_boot, confidence, rng_b)
    boot_widths_ln.append(bw)

# Lognormal CI widths have higher seed-to-seed variance due to skew;
# 1.5× envelope is the minimal bound that absorbs the extra variability
# while still rejecting non-convergent algorithms.
check_true(
    "Log-normal width converges",
    boot_widths_ln[-1] <= boot_widths_ln[0] * 1.5,
)

Heavy-Tailed (t, df=3

heavy = model["heavy_tail"]
rng_data3 = np.random.default_rng(heavy["seed"])
data_ht = rng_data3.standard_t(heavy["df"], size=data_n) * heavy["scale"] + heavy["loc"]

boot_widths_ht = []
for n_boot in replicate_counts:
    rng_b = np.random.default_rng(heavy["seed"] + n_boot)
    bw = bootstrap_ci_width(data_ht, n_boot, confidence, rng_b)
    boot_widths_ht.append(bw)
    if n_boot == replicate_counts[-1]:
        print(f"  Width at n={n_boot}: {bw:.4f} (Gaussian was {boot_widths_g[-1]:.4f})")

check_true(
    "Heavy-tail wider than Gaussian",
    boot_widths_ht[-1] > boot_widths_g[-1] * 0.8,
)

Coverage

n_cov_trials = model["n_trials_coverage"]
true_mean_g = gauss["mu"]

boot_cov = coverage_at_n(
    lambda r: r.normal(gauss["mu"], gauss["sigma"], data_n),
    true_mean_g, bootstrap_ci_width,
    n_cov_trials, 1000, confidence, gauss["seed"] + 9000,
)
rawr_cov = coverage_at_n(
    lambda r: r.normal(gauss["mu"], gauss["sigma"], data_n),
    true_mean_g, rawr_ci_width,
    n_cov_trials, 1000, confidence, gauss["seed"] + 19000,
)
print(f"  Bootstrap coverage (n=1000, {n_cov_trials} trials): {boot_cov:.3f}")
print(f"  RAWR coverage (n=1000, {n_cov_trials} trials):      {rawr_cov:.3f}")

check_true(
    "Bootstrap coverage ≥ 85%",
    boot_cov >= exp["bootstrap_coverage_min"],
)
check_true(
    "RAWR coverage ≥ 82%",
    rawr_cov >= exp["rawr_coverage_min"],
)

Determinism

det_data = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0])
rng1 = np.random.default_rng(7777)
rng2 = np.random.default_rng(7777)
w1 = bootstrap_ci_width(det_data, 500, 0.95, rng1)
w2 = bootstrap_ci_width(det_data, 500, 0.95, rng2)
check_true("Bootstrap deterministic", w1 == w2)

rng3 = np.random.default_rng(8888)
rng4 = np.random.default_rng(8888)
w3 = rawr_ci_width(det_data, 500, 0.95, rng3)
w4 = rawr_ci_width(det_data, 500, 0.95, rng4)
check_true("RAWR deterministic", w3 == w4)

Key Findings

print(f"\n{'=' * 72}")

Both methods converge by ~2000 replicates for Gaussian data

print("   → for most groundSpring experiments, n=2000 is sufficient")

# Results: {pass_count()}/{total_count()} checks passed
print_summary("Exp 013: Resampling Convergence")

Visualization

# Publication-grade summary chart for Exp 013
fig, ax = plt.subplots(figsize=(8, 4))

p, f_count, t = pass_count(), fail_count(), total_count()
ax.barh(['Pass', 'Fail'], [p, f_count], color=[PASS_COLOR, FAIL_COLOR])
ax.set_xlim(0, max(t * 1.15, 1))
ax.set_title('Exp 013: Resampling Convergence Analysis — Validation Results')
ax.set_xlabel('Check Count')
for i, v in enumerate([p, f_count]):
    if v > 0:
        ax.text(v + 0.3, i, str(v), va='center', fontweight='bold')

plt.tight_layout()
plt.savefig(f'/tmp/groundspring_exp013.png', dpi=150, bbox_inches='tight')
plt.show()
print(f'\nResult: {p}/{t} PASS, {f_count}/{t} FAIL')

Provenance & Summary

Provenance chain: Python baseline → Rust validation → barraCuda GPU → metalForge cross-substrate → primal IPC composition

See primals.eco for rendered lab notebooks.