import numpy as np
rng=np.random.default_rng(41)
d=np.arange(180)
# Synthetische dagreeks vanaf01-10-2025, klok UTC+1.
temperatuur=8+7*np.cos(2*np.pi*d/100)+rng.normal(0,1.5,len(d))
kou=np.maximum(0,9.34-temperatuur)
donker=np.clip(8-0.2*temperatuur+rng.normal(0,1,len(d)),3,12)
energie=8+0.38*kou+0.20*donker+rng.normal(0,0.25,len(d))
X=np.column_stack((np.ones(len(d)),kou,donker))
train=slice(0,110); val=slice(110,145); test=slice(145,180)
gem=X[train,1:].mean(axis=0); sd=X[train,1:].std(axis=0)
Z=np.column_stack((np.ones(len(d)),(X[:,1:]-gem)/sd))
def ridge(Z,y,lam):
    # Intercept wordt niet bestraft; solve in plaats van inverse.
    straf=np.diag([0.,lam,lam])
    return np.linalg.solve(Z.T@Z+straf,Z.T@y)
def rmse(y,p): return np.sqrt(np.mean((y-p)**2))
kandidaten=[]
for lam in (0.,0.1,1.,10.,100.):
    b=ridge(Z[train],energie[train],lam)
    kandidaten.append((rmse(energie[val],Z[val]@b),lam,b))
fout,lam,b=min(kandidaten,key=lambda r:r[0])
p=Z@b
print("lambda",lam,"validatieRMSE",round(fout,4))
print("testRMSE",round(rmse(energie[test],p[test]),4),"kWh/dag")
print("basislijn testRMSE",round(rmse(energie[test],np.full(35,energie[train].mean())),4))
fysiek=b[1:]/sd
intercept=b[0]-gem@fysiek
print("coefficienten",round(intercept,4),*[round(x,4) for x in fysiek])
# Kleine rolling-origin CV; schaal alleen op elke trainingsfold.
for eind in (60,80,100):
    a=X[:eind]; y=energie[:eind]
    coef=np.linalg.lstsq(a,y,rcond=None)[0]
    print("fold",eind,"RMSE",round(rmse(energie[eind:eind+20],X[eind:eind+20]@coef),4))
assert rmse(energie[test],p[test])<0.5