sim = {
const mulberry32 = seed => {
let a = seed >>> 0;
return () => {
a |= 0;
a = (a + 0x6D2B79F5) | 0;
let t = Math.imul(a ^ (a >>> 15), 1 | a);
t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
};
};
const normal = rng => {
let u = 0, v = 0;
while (u === 0) u = rng();
while (v === 0) v = rng();
return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);
};
const mean = a => a.reduce((s, x) => s + x, 0) / a.length;
const variance = a => {
const m = mean(a);
return a.reduce((s, x) => s + (x - m) ** 2, 0) / Math.max(1, a.length - 1);
};
const covariance = (x, y) => {
const mx = mean(x), my = mean(y);
return x.reduce((s, xi, i) => s + (xi - mx) * (y[i] - my), 0) / Math.max(1, x.length - 1);
};
const ols = (x, y) => {
const mx = mean(x), my = mean(y);
const sxx = x.reduce((s, xi) => s + (xi - mx) ** 2, 0);
const sxy = x.reduce((s, xi, i) => s + (xi - mx) * (y[i] - my), 0);
const slope = sxx === 0 ? 0 : sxy / sxx;
const intercept = my - slope * mx;
const fitted = x.map(xi => intercept + slope * xi);
const sse = y.reduce((s, yi, i) => s + (yi - fitted[i]) ** 2, 0);
const sst = y.reduce((s, yi) => s + (yi - my) ** 2, 0);
const r2 = sst === 0 ? 0 : 1 - sse / sst;
const rmse = Math.sqrt(sse / y.length);
return {intercept, slope, fitted, r2, rmse, sse};
};
const quantile = (a, q) => {
const s = [...a].sort((x, y) => x - y);
const pos = (s.length - 1) * q;
const lo = Math.floor(pos), hi = Math.ceil(pos);
if (lo === hi) return s[lo];
return s[lo] + (s[hi] - s[lo]) * (pos - lo);
};
const sigmoid = z => 1 / (1 + Math.exp(-z));
const solve2 = (a, b, c, d, e, f) => {
const det = a * d - b * c;
if (Math.abs(det) < 1e-10) return [0, 0];
return [(e * d - b * f) / det, (a * f - e * c) / det];
};
const logisticFit = (x, y) => {
let b0 = 0, b1 = 0;
for (let it = 0; it < 30; it++) {
let g0 = 0, g1 = 0, h00 = 1e-6, h01 = 0, h11 = 1e-6;
for (let i = 0; i < x.length; i++) {
const p = sigmoid(Math.max(-20, Math.min(20, b0 + b1 * x[i])));
const w = Math.max(1e-6, p * (1 - p));
const r = y[i] - p;
g0 += r;
g1 += r * x[i];
h00 += w;
h01 += w * x[i];
h11 += w * x[i] * x[i];
}
const [d0, d1] = solve2(h00, h01, h01, h11, g0, g1);
b0 += d0; b1 += d1;
if (Math.max(Math.abs(d0), Math.abs(d1)) < 1e-7) break;
}
return [b0, b1];
};
const poissonFit = (x, y) => {
let b0 = Math.log(Math.max(0.1, mean(y))), b1 = 0;
for (let it = 0; it < 30; it++) {
let g0 = 0, g1 = 0, h00 = 1e-6, h01 = 0, h11 = 1e-6;
for (let i = 0; i < x.length; i++) {
const eta = Math.max(-8, Math.min(8, b0 + b1 * x[i]));
const mu = Math.exp(eta);
const r = y[i] - mu;
g0 += r;
g1 += r * x[i];
h00 += mu;
h01 += mu * x[i];
h11 += mu * x[i] * x[i];
}
const [d0, d1] = solve2(h00, h01, h01, h11, g0, g1);
b0 += d0; b1 += d1;
if (Math.max(Math.abs(d0), Math.abs(d1)) < 1e-7) break;
}
return [b0, b1];
};
const poissonQuantile = (u, lambda) => {
const lam = Math.max(0.001, Math.min(150, lambda));
let k = 0;
let p = Math.exp(-lam);
let cdf = p;
while (u > cdf && k < 400) {
k += 1;
p *= lam / k;
cdf += p;
}
return k;
};
const binomialPmf = (k, n, p) => {
if (k < 0 || k > n || !Number.isInteger(k)) return 0;
const pp = Math.max(1e-12, Math.min(1 - 1e-12, p));
return Math.exp(
logGamma(n + 1) - logGamma(k + 1) - logGamma(n - k + 1) +
k * Math.log(pp) + (n - k) * Math.log1p(-pp)
);
};
const poissonPmf = (k, lambda) => {
if (k < 0 || !Number.isInteger(k) || lambda <= 0) return 0;
return Math.exp(k * Math.log(lambda) - lambda - logGamma(k + 1));
};
const logGamma = z => {
const c = [
0.9999999999998099, 676.5203681218851, -1259.1392167224028,
771.3234287776531, -176.6150291621406, 12.50734327868691,
-0.1385710952657201, 9.984369578019572e-6, 1.505632735149312e-7
];
if (z < 0.5) return Math.log(Math.PI) - Math.log(Math.sin(Math.PI * z)) - logGamma(1 - z);
let x = c[0];
z -= 1;
for (let i = 1; i < c.length; i++) x += c[i] / (z + i);
const t = z + 7.5;
return 0.5 * Math.log(2 * Math.PI) + (z + 0.5) * Math.log(t) - t + Math.log(x);
};
const betaCF = (a, b, x) => {
const maxIt = 200, eps = 3e-10, tiny = 1e-30;
const qab = a + b, qap = a + 1, qam = a - 1;
let c = 1, d = 1 - qab * x / qap;
if (Math.abs(d) < tiny) d = tiny;
d = 1 / d;
let h = d;
for (let m = 1; m <= maxIt; m++) {
const m2 = 2 * m;
let aa = m * (b - m) * x / ((qam + m2) * (a + m2));
d = 1 + aa * d; if (Math.abs(d) < tiny) d = tiny;
c = 1 + aa / c; if (Math.abs(c) < tiny) c = tiny;
d = 1 / d; h *= d * c;
aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
d = 1 + aa * d; if (Math.abs(d) < tiny) d = tiny;
c = 1 + aa / c; if (Math.abs(c) < tiny) c = tiny;
d = 1 / d;
const del = d * c;
h *= del;
if (Math.abs(del - 1) < eps) break;
}
return h;
};
const regBeta = (x, a, b) => {
if (x <= 0) return 0;
if (x >= 1) return 1;
const bt = Math.exp(logGamma(a + b) - logGamma(a) - logGamma(b) + a * Math.log(x) + b * Math.log1p(-x));
return x < (a + 1) / (a + b + 2)
? bt * betaCF(a, b, x) / a
: 1 - bt * betaCF(b, a, 1 - x) / b;
};
const tTwoSidedP = (t, df) => {
if (!Number.isFinite(t) || df <= 0) return NaN;
return Math.min(1, Math.max(0, regBeta(df / (df + t * t), df / 2, 0.5)));
};
const fUpperP = (f, df1, df2) => {
if (!Number.isFinite(f) || f < 0 || df1 <= 0 || df2 <= 0) return NaN;
return Math.min(1, Math.max(0, regBeta(df2 / (df2 + df1 * f), df2 / 2, df1 / 2)));
};
const erf = x => {
const sign = x < 0 ? -1 : 1;
x = Math.abs(x);
const a1=0.254829592,a2=-0.284496736,a3=1.421413741,a4=-1.453152027,a5=1.061405429,p=0.3275911;
const t = 1 / (1 + p * x);
const y = 1 - (((((a5*t+a4)*t+a3)*t+a2)*t+a1)*t) * Math.exp(-x*x);
return sign * y;
};
const zTwoSidedP = z => 2 * (1 - 0.5 * (1 + erf(Math.abs(z) / Math.sqrt(2))));
const slopeSEFromWeights = (x, w) => {
const h00 = w.reduce((s, wi) => s + wi, 0);
const h01 = w.reduce((s, wi, i) => s + wi * x[i], 0);
const h11 = w.reduce((s, wi, i) => s + wi * x[i] * x[i], 0);
const det = h00 * h11 - h01 * h01;
return det <= 0 ? NaN : Math.sqrt(h00 / det);
};
const fmtP = p => !Number.isFinite(p) ? "not available" : (p < 0.001 ? "< .001" : `= ${p.toFixed(3)}`);
const olsMatrix = (X, y) => {
const n = X.length, p = X[0].length;
const xtx = Array.from({length:p},()=>Array(p).fill(0));
const xty = Array(p).fill(0);
for (let i=0;i<n;i++) {
for (let j=0;j<p;j++) {
xty[j] += X[i][j] * y[i];
for (let k=0;k<p;k++) xtx[j][k] += X[i][j] * X[i][k];
}
}
const aug = xtx.map((row,i)=>[...row,...Array.from({length:p},(_,j)=>i===j?1:0)]);
for (let col=0;col<p;col++) {
let pivot=col;
for (let r=col+1;r<p;r++) if (Math.abs(aug[r][col])>Math.abs(aug[pivot][col])) pivot=r;
if (Math.abs(aug[pivot][col])<1e-12) return null;
[aug[col],aug[pivot]]=[aug[pivot],aug[col]];
const div=aug[col][col];
for (let j=0;j<2*p;j++) aug[col][j]/=div;
for (let r=0;r<p;r++) if (r!==col) {
const f=aug[r][col];
for (let j=0;j<2*p;j++) aug[r][j]-=f*aug[col][j];
}
}
const inv=aug.map(row=>row.slice(p));
const beta=inv.map(row=>row.reduce((s,v,j)=>s+v*xty[j],0));
const fitted=X.map(row=>row.reduce((s,v,j)=>s+v*beta[j],0));
const sse=y.reduce((s,yi,i)=>s+(yi-fitted[i])**2,0);
const df=n-p, mse=sse/df;
const se=inv.map((row,j)=>Math.sqrt(Math.max(0,mse*row[j])));
const t=beta.map((b,j)=>b/se[j]);
const pValues=t.map(v=>tTwoSidedP(v,df));
return {beta,se,t,p:pValues,fitted,sse,df};
};
return {
mulberry32, normal, mean, variance, covariance, ols, quantile,
sigmoid, logisticFit, poissonFit, poissonQuantile,
binomialPmf, poissonPmf,
tTwoSidedP, fUpperP, zTwoSidedP, slopeSEFromWeights, fmtP, olsMatrix
};
}
