# 10,000 simulated studies with the rhythm.csv design (12 time points, ZT0 to ZT44, 3 mice per time point and genotype,
# residual SD 0.4). Python 3.13, NumPy, SciPy. Seed fixed before running.
import numpy as np
from scipy import stats
rng = np.random.default_rng(20261111); R = 10_000; w = 2 * np.pi / 24
zt = np.repeat(np.arange(0, 48, 4), 3)
def sse(X, y): b, *_ = np.linalg.lstsq(X, y, rcond=None); r = y - X @ b; return r @ r
def study(aw, ak, pw, pk):
    yw = 5 + aw * np.cos(w * (zt - pw)) + rng.normal(0, 0.4, zt.size)
    yk = 5 + ak * np.cos(w * (zt - pk)) + rng.normal(0, 0.4, zt.size)
    y = np.r_[yw, yk]; t = np.r_[zt, zt]; k = np.r_[np.zeros(zt.size), np.ones(zt.size)]
    c, s = np.cos(w * t), np.sin(w * t); one = np.ones_like(t)
    full = np.column_stack([one, c, s, k, k * c, k * s]); red = np.column_stack([one, c, s, k])
    sf, sr = sse(full, y), sse(red, y); dfe = y.size - 6
    p_f = stats.f.sf(((sr - sf) / 2) / (sf / dfe), 2, dfe)
    p_t = min(stats.ttest_ind(yk[zt == z], yw[zt == z], equal_var=True).pvalue for z in np.unique(zt))
    return p_f < 0.05, p_t < 0.05, p_t < 0.05 / 12          # Bonferroni: any test below 0.05/12
for lab, args in (("Same rhythm in both genotypes", (1.0, 1.0, 10, 10)), ("KO amplitude halved (1.0 vs 0.5)", (1.0, 0.5, 10, 10)),
                  ("KO peak 3 h later", (1.0, 1.0, 10, 13))):
    r = np.array([study(*args) for _ in range(R)])
    print(f"{lab}: cosinor F-test for a difference p<0.05 {r[:, 0].mean():.1%}; at least one of 12 time-point t-tests p<0.05 {r[:, 1].mean():.1%}; with Bonferroni correction {r[:, 2].mean():.1%}")