1import json
 2import numpy as np
 3from spectral_experiment import recovery_sweep
 4
 5# For y=rho*(z^2-1)+eps, E[(y-Ey)(z^2-1)] = 2*rho exactly,
 6# since E[(z^2-1)^2]=2 and eps is independent.  This is the rank-one
 7# population spike predicted by Stein's identity.
 8def population_scaling(rhos=(.125,.25,.5,1.,2.), N=500000, noise=.35):
 9    rng=np.random.default_rng(9917)
10    z=rng.normal(size=N); eps=rng.normal(size=N)
11    rows=[]
12    for rho in rhos:
13        y=rho*(z*z-1)+noise*eps
14        # No clipping here: this isolates the population identity.
15        alpha=np.mean((y-y.mean())*(z*z-1))
16        rows.append({'rho':rho,'predicted_2rho':2*rho,'measured_alpha':float(alpha),
17                     'relative_error':float(abs(alpha-2*rho)/(2*rho))})
18    # Independent labels imply zero signal coefficient.
19    y0=rng.normal(size=N)
20    null=np.mean((y0-y0.mean())*(z*z-1))
21    return rows,float(null)
22
23def transition_prediction(d=96, rho=.5, noise=.35):
24    # Leading-order signed-Wishart heuristic: alpha=2rho, weighted-noise
25    # variance s2=E[(y-Ey)^2]. The predicted outlier boundary is
26    # gamma=d/n < alpha^2/s2, i.e. n/d > s2/alpha^2.
27    alpha=2*rho
28    s2=2*rho*rho+noise*noise
29    critical=s2/(alpha*alpha)
30    return {'rho':rho,'alpha_predicted':alpha,'weight_variance_predicted':s2,
31            'predicted_critical_n_over_d':critical}
32
33def main():
34    scaling,null=population_scaling()
35    trans=[]
36    for rho in (.5,1.): trans.append(transition_prediction(rho=rho))
37    # Reuse the actual estimator sweep and identify first ratio with overlap
38    # clearly above the null 1/d baseline (a conservative finite-size onset).
39    rec=recovery_sweep(d=96,rhos=(.5,1.),ratios=(.5,1,2,4,8),trials=8)
40    onset=[]
41    for rho in (.5,1.):
42        rr=[x for x in rec if x['rho']==rho]
43        hit=next((x['n_over_d'] for x in rr if x['overlap_mean']>4/96),None)
44        onset.append({'rho':rho,'observed_onset_n_over_d':hit,
45                      'null_overlap_prediction':1/96})
46    out={'population_scaling':scaling,'null_alpha':null,
47         'transition_predictions':trans,'observed_recovery_onsets':onset,
48         'recovery_sweep':rec}
49    with open('mechanism_results.json','w') as f: json.dump(out,f,indent=2)
50    print(json.dumps(out,indent=2))
51if __name__=='__main__': main()