# 10,000 simulated experiments as in make_data.py (treatment lowers body weight by 4 g; no direct effect on the
# heart; heart = 5.2 * BW^0.75). Share of runs in which each analysis declares a heart effect (p < 0.05).
import numpy as np
from scipy import stats
rng = np.random.default_rng(11)
runs, n, hit_ratio, hit_ancova = 10_000, 10, 0, 0
g = np.r_[np.zeros(n), np.ones(n)]
for _ in range(runs):
    bw = np.r_[28 + rng.normal(0, 2, n), 24 + rng.normal(0, 2, n)]
    heart = 5.2 * bw ** 0.75 * np.exp(rng.normal(0, 0.06, 2 * n))
    r = heart / bw
    hit_ratio += stats.ttest_ind(r[g == 1], r[g == 0], equal_var=False).pvalue < 0.05
    X = np.c_[np.ones(2 * n), g, np.log(bw)]
    beta, res, *_ = np.linalg.lstsq(X, np.log(heart), rcond=None)
    sigma2 = res[0] / (2 * n - 3)
    se = np.sqrt(sigma2 * np.linalg.inv(X.T @ X)[1, 1])
    hit_ancova += 2 * stats.t.sf(abs(beta[1] / se), 2 * n - 3) < 0.05
print(f"runs = {runs}; heart effect declared: heart/body-weight ratio {hit_ratio / runs:.1%}, ANCOVA on logs {hit_ancova / runs:.1%}")