"""Authored OLS fixture, known sigma assumption: mean confidence versus new observation normal interval."""
from fractions import Fraction as F
from statistics import NormalDist
import math,json
X=[F(-2),F(-1),F(0),F(1),F(2)];Y=[x+1 for x in X]
def demo():
    n=len(X);xm=sum(X,F(0))/n;ym=sum(Y,F(0))/n
    sxx=sum(((x-xm)**2 for x in X),F(0));sxy=sum(((x-xm)*(y-ym) for x,y in zip(X,Y)),F(0))
    slope=sxy/sxx;intercept=ym-slope*xm
    assert slope==intercept==1 and sxx==10
    sigma2=F(1);z=NormalDist().inv_cdf(.975);rows=[]
    for x in [F(0),F(3)]:
        mean=intercept+slope*x;h=F(1,n)+(x-xm)**2/sxx
        vm=sigma2*h;vo=sigma2*(h+1)
        half_mean=z*math.sqrt(float(vm));half_observation=z*math.sqrt(float(vo))
        assert vo-vm==sigma2 and half_observation>half_mean
        rows.append({'x':str(x),'predicted_mean':str(mean),'leverage_h':str(h),'mean_variance':str(vm),'new_observation_variance':str(vo),'mean_confidence_interval':[float(mean)-half_mean,float(mean)+half_mean],'new_observation_prediction_interval':[float(mean)-half_observation,float(mean)+half_observation],'outside_training_x_range':x<min(X) or x>max(X)})
    residual_sse=sum(((y-intercept-slope*x)**2 for x,y in zip(X,Y)),F(0));assert residual_sse==0
    return {'authored_x':[str(x) for x in X],'authored_y':[str(y) for y in Y],'actual_ols_intercept':str(intercept),'actual_ols_slope':str(slope),'Sxx':str(sxx),'known_noise_variance_assumed_not_estimated':str(sigma2),'authored_fixture_residual_sse':str(residual_sse),'model_assumptions':['correct linear mean','independent Gaussian errors conditional on fixed X','known sigma=1, not estimated from authored residuals','new error independent of fitted mean'],'interval_level':.95,'critical_distribution':'standard normal because sigma assumed known; not estimated-variance t','normal_quantile':z,'rows':rows,'empirical_coverage_measured':False,'statsmodels_executed':False}
if __name__=='__main__':print(json.dumps(demo(),ensure_ascii=False,indent=2))