Marchenko-Pastur: Eigenvalues vs Theory

const N = 80, T = 200, Q = N / T;
const LM = (1 - Math.sqrt(Q)) ** 2, LP = (1 + Math.sqrt(Q)) ** 2;
const NBINS = 50;
const BIN_LO = Math.max(0, LM - 0.15), BIN_HI = LP + 0.15;
const BW = (BIN_HI - BIN_LO) / NBINS;
const hist = new Float64Array(NBINS);
const C = []; for (let i = 0; i < N; i++) C.push(new Float64Array(N));
const X = []; for (let i = 0; i < N; i++) X.push(new Float64Array(T));
let smooth = new Float64Array(NBINS);

function nrand() {
  let s = 0;
  for (let k = 0; k < 12; k++) s += Math.random();
  return s - 6;
}

function sampleCov() {
  for (let i = 0; i < N; i++) {
    const Xi = X[i];
    for (let t = 0; t < T; t++) Xi[t] = nrand();
  }
  for (let i = 0; i < N; i++) {
    const Xi = X[i], Ci = C[i];
    for (let j = i; j < N; j++) {
      const Xj = X[j];
      let s = 0;
      for (let t = 0; t < T; t++) s += Xi[t] * Xj[t];
      s /= T;
      Ci[j] = s;
      C[j][i] = s;
    }
  }
}

function jacobiEigs(A, sweeps) {
  const n = A.length;
  const M = []; for (let i = 0; i < n; i++) M.push(Float64Array.from(A[i]));
  for (let sw = 0; sw < sweeps; sw++) {
    let off = 0;
    for (let p = 0; p < n - 1; p++) {
      for (let q = p + 1; q < n; q++) {
        const apq = M[p][q];
        if (Math.abs(apq) < 1e-10) continue;
        off += apq * apq;
        const app = M[p][p], aqq = M[q][q];
        const theta = (aqq - app) / (2 * apq);
        const t = (theta >= 0 ? 1 : -1) / (Math.abs(theta) + Math.sqrt(1 + theta * theta));
        const c = 1 / Math.sqrt(1 + t * t), s = t * c;
        M[p][p] = app - t * apq;
        M[q][q] = aqq + t * apq;
        M[p][q] = 0; M[q][p] = 0;
        for (let r = 0; r < n; r++) {
          if (r !== p && r !== q) {
            const arp = M[r][p], arq = M[r][q];
            M[r][p] = c * arp - s * arq;
            M[r][q] = s * arp + c * arq;
            M[p][r] = M[r][p];
            M[q][r] = M[r][q];
          }
        }
      }
    }
    if (off < 1e-10) break;
  }
  const eig = new Float64Array(n);
  for (let i = 0; i < n; i++) eig[i] = M[i][i];
  return eig;
}

function mpDensity(l) {
  if (l <= LM || l >= LP) return 0;
  return Math.sqrt((LP - l) * (l - LM)) / (2 * Math.PI * Q * l);
}

function init({ ctx, width, height }) {
  ctx.fillStyle = '#0a0a12';
  ctx.fillRect(0, 0, width, height);
}

function tick({ ctx, width, height, frame }) {
  sampleCov();
  const eig = jacobiEigs(C, 40);
  hist.fill(0);
  for (let i = 0; i < N; i++) {
    const b = Math.floor((eig[i] - BIN_LO) / BW);
    if (b >= 0 && b < NBINS) hist[b] += 1;
  }
  const norm = 1 / (N * BW);
  const a = frame < 3 ? 1 : 0.18;
  for (let i = 0; i < NBINS; i++) smooth[i] = smooth[i] * (1 - a) + hist[i] * norm * a;
  ctx.fillStyle = 'rgba(10,10,18,1)';
  ctx.fillRect(0, 0, width, height);
  const padL = 60, padR = 30, padT = 50, padB = 50;
  const W = width - padL - padR, H = height - padT - padB;
  let yMax = 0;
  for (let i = 0; i < NBINS; i++) if (smooth[i] > yMax) yMax = smooth[i];
  const mpPeak = mpDensity((LM + LP) / 2);
  yMax = Math.max(yMax, mpPeak) * 1.2;
  const xOf = l => padL + ((l - BIN_LO) / (BIN_HI - BIN_LO)) * W;
  const yOf = v => padT + H - (v / yMax) * H;
  ctx.strokeStyle = '#23243a';
  ctx.lineWidth = 1;
  ctx.beginPath();
  for (let g = 0; g <= 5; g++) {
    const y = padT + (H * g) / 5;
    ctx.moveTo(padL, y); ctx.lineTo(padL + W, y);
  }
  ctx.stroke();
  for (let i = 0; i < NBINS; i++) {
    const x0 = xOf(BIN_LO + i * BW), x1 = xOf(BIN_LO + (i + 1) * BW);
    const y0 = yOf(smooth[i]);
    const h = padT + H - y0;
    const grad = ctx.createLinearGradient(0, y0, 0, padT + H);
    grad.addColorStop(0, '#5ad1ff');
    grad.addColorStop(1, '#1a4d7a');
    ctx.fillStyle = grad;
    ctx.fillRect(x0 + 1, y0, x1 - x0 - 2, h);
  }
  ctx.strokeStyle = '#ff7a3d';
  ctx.lineWidth = 2.5;
  ctx.beginPath();
  let started = false;
  const STEPS = 240;
  for (let k = 0; k <= STEPS; k++) {
    const l = BIN_LO + (BIN_HI - BIN_LO) * (k / STEPS);
    const rho = mpDensity(l);
    if (rho > 0) {
      const x = xOf(l), y = yOf(rho);
      if (!started) { ctx.moveTo(x, y); started = true; } else ctx.lineTo(x, y);
    } else started = false;
  }
  ctx.stroke();
  ctx.strokeStyle = 'rgba(255,122,61,0.45)';
  ctx.setLineDash([4, 4]);
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(xOf(LM), padT); ctx.lineTo(xOf(LM), padT + H);
  ctx.moveTo(xOf(LP), padT); ctx.lineTo(xOf(LP), padT + H);
  ctx.stroke();
  ctx.setLineDash([]);
  ctx.fillStyle = '#e8ecf4';
  ctx.font = '14px ui-sans-serif, system-ui';
  ctx.fillText('Marchenko-Pastur: empirical vs theory', padL, 26);
  ctx.fillStyle = '#9aa3b8';
  ctx.font = '11px ui-monospace, monospace';
  ctx.fillText(`N=${N}  T=${T}  q=${Q.toFixed(2)}`, padL, 42);
  ctx.fillText(`lambda- = ${LM.toFixed(3)}`, xOf(LM) + 4, padT + 12);
  ctx.fillText(`lambda+ = ${LP.toFixed(3)}`, xOf(LP) - 78, padT + 12);
  ctx.fillStyle = '#9aa3b8';
  ctx.font = '10px ui-monospace, monospace';
  for (let g = 0; g <= 5; g++) {
    const l = BIN_LO + ((BIN_HI - BIN_LO) * g) / 5;
    const x = xOf(l);
    ctx.fillText(l.toFixed(2), x - 10, padT + H + 14);
  }
}