"""Case 007: same planning assumptions and hypothetical summaries as R.
Requires NumPy and SciPy. n is per group, not total sample size.
"""
import numpy as np
import scipy
from scipy.stats import t, nct
from scipy.optimize import brentq
sd_assumed=10.0
important_difference=3.0
alpha=.05
def design_power(n,delta):
    df=2*n-2
    critical=t.ppf(1-alpha/2,df)
    noncentrality=delta/(sd_assumed*np.sqrt(2/n))
    return nct.sf(critical,df,noncentrality)+nct.sf(critical,df,-noncentrality)
def summarize(n):
    se=sd_assumed*np.sqrt(2/n)
    half_width=t.ppf(1-alpha/2,2*n-2)*se
    # Assumes sample SDs = 10 and observed mean difference = 3.
    return np.array([n,2*n,se,3-half_width,3+half_width,design_power(n,important_difference)])
planned,achieved=summarize(180),summarize(60)
print("per_group total SE CI_low CI_high power_for_prespecified_difference")
print("planned:",np.round(planned,6))
print("achieved:",np.round(achieved,6))
mde=brentq(lambda delta:design_power(60,delta)-.80,0,10)
print("difference_for_80_percent_power:",round(mde,6))
print("SE_ratio:",round(achieved[2]/planned[2],6))
assert np.isclose(achieved[2]/planned[2],np.sqrt(3),atol=1e-10)
assert planned[3]>0 and achieved[3]<0
assert achieved[5]<planned[5] and mde>important_difference
print("All assertions passed; NumPy",np.__version__,"SciPy",scipy.__version__)