"""Larger check of Welch 95% CI coverage at n = 3 per group (normal data, equal SDs).
The 10,000-run figures in example.py (Python) and example.R differ by about 0.7 points; this
200,000-run vectorised check settles the value. Python 3.13, NumPy, SciPy 1.18.1.
"""
import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
N, n = 200_000, 3
x = rng.normal(0, 1, (N, n))
y = rng.normal(1, 1, (N, n))
v1 = x.var(axis=1, ddof=1) / n
v2 = y.var(axis=1, ddof=1) / n
se = np.sqrt(v1 + v2)
df = (v1 + v2) ** 2 / (v1**2 / (n - 1) + v2**2 / (n - 1))
d = y.mean(axis=1) - x.mean(axis=1)
q = stats.t.ppf(0.975, df)
print("coverage n=3, 200000 runs:", round(float(np.mean((d - q * se <= 1) & (1 <= d + q * se))), 4))