1import json
  2import numpy as np
  3from pathlib import Path
  4
  5EPS = 1e-10
  6
  7def corr_eigh(X):
  8    X = X - X.mean(0, keepdims=True)
  9    s = X.std(0, ddof=1) + 1e-10
 10    C = (X / s).T @ (X / s) / X.shape[0]
 11    w, U = np.linalg.eigh(C)
 12    order = np.argsort(w)[::-1]
 13    return w[order], U[:, order]
 14
 15def band_statistics(X, bands):
 16    w, U = corr_eigh(X)
 17    Z = (X - X.mean(0)) @ U
 18    vals, us = [], []
 19    for a, b in bands:
 20        z = Z[:, a:b]
 21        v = np.mean(z*z, axis=0).mean()
 22        # aggregate fourth cumulant over modes in the band
 23        k4 = (np.mean(z**4, axis=0) - 3*np.mean(z*z, axis=0)**2).mean()
 24        vals.append(np.mean(w[a:b]))
 25        us.append(-k4/(v*v + EPS))
 26    vals, us = np.asarray(vals), np.asarray(us)
 27    d = np.diff(np.log(np.abs(us)+1e-7)) / np.diff(np.log(vals+1e-7))
 28    return w, U, vals, us, d
 29
 30def rg_mask(w, us, bands, retain):
 31    # Practical RG proxy: a persistent change in log quartic coupling identifies
 32    # the signal side; rank bands by slope magnitude, then keep whole bands.
 33    scale = np.asarray([w[a:b].mean() for a,b in bands])
 34    d = np.diff(np.log(np.abs(us)+1e-7)) / np.diff(np.log(scale+1e-7))
 35    score = np.r_[np.abs(d), np.abs(d[-1]) if len(d) else 0]
 36    # retain contiguous/whole bands with the strongest scale dependence
 37    chosen = np.argsort(score)[::-1][:max(1, int(np.ceil(retain / (bands[0][1]-bands[0][0]))))]
 38    mask = np.zeros(len(w), dtype=bool)
 39    for r in chosen:
 40        a,b = bands[r]; mask[a:b] = True
 41    # exact width for fair comparison
 42    ids = np.flatnonzero(mask)
 43    if len(ids) > retain: mask[ids[retain:]] = False
 44    return mask, d, score
 45
 46def make_controlled(n=12000, p=40, seed=0, s=1.5, signal_bands=(4,7), strength=1.0):
 47    rng=np.random.default_rng(seed)
 48    lam=np.linspace(2.0,0.25,p)
 49    X=rng.normal(size=(n,p))*np.sqrt(lam)
 50    # Non-Gaussian modes distributed in the middle/bulk. Laplace has excess kurtosis 3.
 51    a,b=signal_bands
 52    for j in range(a,b):
 53        z=rng.laplace(size=n)/np.sqrt(2) # variance 1, excess kurtosis 3
 54        # preserve covariance eigenvalue and impose u4 approximately -strength*lambda^s
 55        X[:,j]=np.sqrt(lam[j])*((1-strength)*rng.normal(size=n)+strength*z)
 56    return X, lam
 57
 58def verification():
 59    # Prediction 1: the Gaussian excess-kurtosis estimator has zero mean and
 60    # standard deviation proportional to N^-1/2.
 61    ns=[500,2000,8000]
 62    noise=[]
 63    for n in ns:
 64        vals=[]
 65        for seed in range(30):
 66            x=np.random.default_rng(seed).normal(size=(n,8))
 67            k=np.mean(x**4,0)-3*np.mean(x*x,0)**2
 68            vals.append(np.mean(k))
 69        noise.append(float(np.std(vals)))
 70    scaling=float(np.polyfit(np.log(ns),np.log(noise),1)[0])
 71
 72    # Prediction 2: if excess kurtosis kappa4 = c*lambda^(2+s), then
 73    # u4=-kappa4/variance^2 has logarithmic slope s.  Use known independent
 74    # modes so eigendecomposition noise cannot obscure this algebraic check.
 75    slopes=[]
 76    rng=np.random.default_rng(4)
 77    lam=np.linspace(0.5,2.0,40)
 78    for s in [0.5,1.0,2.0]:
 79        us=[]
 80        for l in lam:
 81            alpha=0.18*l**s/(2.0**s) # mixture excess kurtosis = 3 alpha
 82            choose=rng.random(12000)<alpha
 83            z=rng.normal(size=12000)
 84            lap=rng.laplace(size=12000)/np.sqrt(2)
 85            z=np.where(choose,lap,z)*np.sqrt(l)
 86            v=np.mean(z*z); k=np.mean(z**4)-3*v*v
 87            us.append(abs(-k/(v*v+EPS)))
 88        slopes.append(float(np.polyfit(np.log(lam),np.log(np.asarray(us)+1e-7),1)[0]))
 89
 90    # Prediction 3: for a fixed variance mixture, excess kurtosis (and |u4|)
 91    # grows approximately linearly with the non-Gaussian mixture strength.
 92    strengths=np.linspace(0,1,5)
 93    observed=[]
 94    for st in strengths:
 95        vals=[]
 96        for seed in range(10):
 97            rr=np.random.default_rng(seed); choose=rr.random(10000)<(0.22*st)
 98            z=np.where(choose,rr.laplace(size=10000)/np.sqrt(2),rr.normal(size=10000))
 99            v=np.mean(z*z); vals.append(abs((np.mean(z**4)-3*v*v)/(v*v+EPS)))
100        observed.append(float(np.mean(vals)))
101    linear=float(np.polyfit(strengths,observed,1)[0])
102    return {'noise_std_by_N':dict(zip(ns,noise)), 'noise_loglog_slope_observed':scaling,
103            'noise_prediction':-0.5, 'controlled_slope_prediction':[.5,1.,2.],
104            'controlled_slope_observed':slopes, 'strengths':strengths.tolist(),
105            'abs_u4_by_strength':observed, 'u4_strength_linear_slope':linear,
106            'u4_strength_prediction':'linear'}
107
108def classification(seed=20):
109    rng=np.random.default_rng(seed); p=40; ntr=3000; nte=1500
110    # Signal is extensive-rank: classes differ in several middle eigenmodes,
111    # while leading modes are high-variance nuisance directions.
112    lam=np.linspace(2.5,.3,p); sig=np.arange(12,28)
113    def gen(n):
114        y=rng.integers(0,2,n); X=rng.normal(size=(n,p))*np.sqrt(lam)
115        X[:,sig] += (2*y[:,None]-1)*0.42*np.sqrt(lam[sig])[None,:]
116        return X,y
117    cal,_=gen(12000); Xtr,ytr=gen(ntr); Xte,yte=gen(nte)
118    bands=[(i,i+4) for i in range(0,p,4)]
119    w,U,bv,u,d=band_statistics(cal,bands)
120    retain=16
121    rmask,_,score=rg_mask(w,u,bands,retain)
122    top=np.zeros(p,bool); top[:retain]=1
123    def acc(mask):
124        # Linear discriminant in retained PCA coordinates, fit without sklearn.
125        ztr=(Xtr-Xtr.mean(0))@U[:,mask]; zte=(Xte-Xtr.mean(0))@U[:,mask]
126        a=ztr[ytr==1].mean(0)-ztr[ytr==0].mean(0)
127        c=.5*(ztr[ytr==1].mean(0)+ztr[ytr==0].mean(0))
128        pred=( (zte-c)@a > 0 )
129        return float(np.mean(pred==yte))
130    return {'retained_width':retain,'rg_selected_bands':np.flatnonzero(rmask.reshape(-1,4).any(1)).tolist(),
131            'top_selected_bands':list(range(4)), 'rg_accuracy':acc(rmask),
132            'top_pca_accuracy':acc(top), 'quartic_u4':u.tolist(), 'slopes':d.tolist()}
133
134def main():
135    out={'verification':verification(),'classification':classification()}
136    Path('results.json').write_text(json.dumps(out,indent=2))
137    print(json.dumps(out,indent=2))
138if __name__=='__main__': main()