# Same deterministic synthetic data and OLS contrasts as R.
import numpy as np
from scipy.stats import t
x=np.tile(np.arange(8.),2); g=np.repeat([0.,1.],8)
X=np.column_stack([np.ones(16),x,g,x*g])
z=np.sin(np.arange(1.,17.))
e=(z-X@np.linalg.lstsq(X,z,rcond=None)[0])*2
y=X@np.array([40.,2.,5.,-1.])+e
b=np.linalg.lstsq(X,y,rcond=None)[0]; df=12
v=(np.sum((y-X@b)**2)/df)*np.linalg.inv(X.T@X)
print(b)
for label,L in [('slope g=0',[0,1,0,0]),('slope g=1',[0,1,0,1]),('slope difference',[0,0,0,1]),('group difference x= 2',[0,0,1,2]),('group difference x= 6',[0,0,1,6])]:
    L=np.array(L); est=L@b; se=np.sqrt(L@v@L); ci=est+np.array([-1,1])*t.ppf(.975,df)*se
    print(f'{label} estimate={est:.6f} SE={se:.6f} CI=[{ci[0]:.6f},{ci[1]:.6f}]')
assert np.allclose(b,[40,2,5,-1])