# DASS Analysis Clinic Case 002; Python dependencies: numpy, scipy.
from pathlib import Path
import numpy as np
from scipy import stats, optimize
HERE = Path(__file__).resolve().parent

def ols(X, y):
    b = np.linalg.lstsq(X, y, rcond=None)[0]
    residual = y - X @ b
    df = len(y) - X.shape[1]
    cov = (residual @ residual / df) * np.linalg.inv(X.T @ X)
    se = np.sqrt(np.diag(cov))
    ci = np.column_stack((b-stats.t.ppf(.975, df)*se,b+stats.t.ppf(.975,df)*se))
    return b, cov, ci

def power(n, d):
    # Matches R power.t.test(strict=FALSE): rejection tail in effect direction.
    return stats.nct.sf(stats.t.ppf(.975,2*n-2),2*n-2,d*np.sqrt(n/2))

p=np.array([.004,.011,.019,.032,.041,.09,.21,.48,.73])
def adjust(p):
 order=np.argsort(p); ordered=p[order]; n=len(p)
 holm=np.minimum(1,np.maximum.accumulate(ordered*(n-np.arange(n))))
 bh=np.minimum(1,np.minimum.accumulate((ordered*n/np.arange(1,n+1))[::-1])[::-1])
 restore=np.argsort(order)
 return np.minimum(1,n*p),holm[restore],bh[restore]
for k in [1,3,5,10,20]: print(f'Independent null tests={k}: FWER={1-.95**k:.3f}')
rng=np.random.default_rng(2026)
print('Python simulation FWER:',np.mean(np.any(rng.random((20000,10))<.05,axis=1)))
bonf,holm,bh=adjust(p)
print('raw / Bonferroni / Holm / BH:'); print(np.column_stack((p,bonf,holm,bh)))
print('Significant counts:',*(np.sum(q<.05) for q in [p,bonf,holm,bh]))
print('Secondary family:'); print(np.column_stack((p[1:],*adjust(p[1:])[1:])))
assert [np.sum(q<.05) for q in [p,bonf,holm,bh]]==[5,1,1,2]