# Exact E[V_CR0]/Var(ybar) under the DGP vs Monte Carlo estimate (z-score with MC SE)
import json, numpy as np, pandas as pd
res=json.load(open("out/results.json")); d=pd.read_csv("data/design.csv")
for _,r in d.iterrows():
    n=np.array([int(x) for x in r.cluster_sizes.split("|")]); G=len(n); N=n.sum(); icc=r.icc
    gid=np.repeat(np.arange(G),n); Z=(gid[:,None]==np.arange(G)).astype(float)
    Sig=icc*Z@Z.T+(1-icc)*np.eye(N); A=Z.T-np.outer(n,np.ones(N))/N
    ev0=np.trace(A@Sig@A.T)/N**2; var=(icc*np.sum(n**2)+(1-icc)*N)/N**2
    k=f"cr0_G{G}_icc{icc:.1f}"; mc=res["rel_bias_var"][k]; se=res["mc_se_rel_bias_var"][k]
    print(k, "exact %.4f  MC %.4f  z=%.2f"%(ev0/var, mc, (mc-ev0/var)/se))