spb/uqo-eval

1// UQO Éval — lois de probabilité pour l'inférence (t, F, χ², normale) sans dépendance.2// Implémentations classiques (Numerical Recipes) : fonction bêta incomplète régularisée et gamma incomplète.34function lnGamma(x: number): number {5  const cof = [76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5];6  let y = x;7  const tmp = x + 5.5 - (x + 0.5) * Math.log(x + 5.5);8  let ser = 1.000000000190015;9  for (let j = 0; j < 6; j++) ser += cof[j] / ++y;10  return -tmp + Math.log((2.5066282746310005 * ser) / x);11}1213function betacf(a: number, b: number, x: number): number {14  const MAXIT = 300;15  const EPS = 3e-14;16  const FPMIN = 1e-300;17  const qab = a + b;18  const qap = a + 1;19  const qam = a - 1;20  let c = 1;21  let d = 1 - (qab * x) / qap;22  if (Math.abs(d) < FPMIN) d = FPMIN;23  d = 1 / d;24  let h = d;25  for (let m = 1; m <= MAXIT; m++) {26    const m2 = 2 * m;27    let aa = (m * (b - m) * x) / ((qam + m2) * (a + m2));28    d = 1 + aa * d;29    if (Math.abs(d) < FPMIN) d = FPMIN;30    c = 1 + aa / c;31    if (Math.abs(c) < FPMIN) c = FPMIN;32    d = 1 / d;33    h *= d * c;34    aa = (-(a + m) * (qab + m) * x) / ((a + m2) * (qap + m2));35    d = 1 + aa * d;36    if (Math.abs(d) < FPMIN) d = FPMIN;37    c = 1 + aa / c;38    if (Math.abs(c) < FPMIN) c = FPMIN;39    d = 1 / d;40    const del = d * c;41    h *= del;42    if (Math.abs(del - 1) < EPS) break;43  }44  return h;45}4647/** I_x(a, b) — bêta incomplète régularisée. */48export function betaInc(a: number, b: number, x: number): number {49  if (x <= 0) return 0;50  if (x >= 1) return 1;51  const bt = Math.exp(lnGamma(a + b) - lnGamma(a) - lnGamma(b) + a * Math.log(x) + b * Math.log(1 - x));52  if (x < (a + 1) / (a + b + 2)) return (bt * betacf(a, b, x)) / a;53  return 1 - (bt * betacf(b, a, 1 - x)) / b;54}5556/** P(a, x) — gamma incomplète régularisée. */57export function gammaInc(a: number, x: number): number {58  if (x <= 0) return 0;59  if (x < a + 1) {60    let ap = a;61    let sum = 1 / a;62    let del = sum;63    for (let n = 0; n < 500; n++) {64      ap += 1;65      del *= x / ap;66      sum += del;67      if (Math.abs(del) < Math.abs(sum) * 3e-14) break;68    }69    return sum * Math.exp(-x + a * Math.log(x) - lnGamma(a));70  }71  // fraction continue (Lentz)72  const FPMIN = 1e-300;73  let b = x + 1 - a;74  let c = 1 / FPMIN;75  let d = 1 / b;76  let h = d;77  for (let i = 1; i < 500; i++) {78    const an = -i * (i - a);79    b += 2;80    d = an * d + b;81    if (Math.abs(d) < FPMIN) d = FPMIN;82    c = b + an / c;83    if (Math.abs(c) < FPMIN) c = FPMIN;84    d = 1 / d;85    const del = d * c;86    h *= del;87    if (Math.abs(del - 1) < 3e-14) break;88  }89  return 1 - Math.exp(-x + a * Math.log(x) - lnGamma(a)) * h;90}9192/** Φ(z). */93export function normalCdf(z: number): number {94  return 0.5 * (1 + erf(z / Math.SQRT2));95}9697function erf(x: number): number {98  // Abramowitz-Stegun 7.1.26 raffiné (erreur < 1.5e-7) — suffisant pour des p-valeurs99  const t = 1 / (1 + 0.3275911 * Math.abs(x));100  const y = 1 - ((((1.061405429 * t - 1.453152027) * t + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t * Math.exp(-x * x);101  return x >= 0 ? y : -y;102}103104/** Quantile de la normale standard (Acklam). */105export function normalQuantile(p: number): number {106  if (p <= 0) return -Infinity;107  if (p >= 1) return Infinity;108  const a = [-3.969683028665376e1, 2.209460984245205e2, -2.759285104469687e2, 1.38357751867269e2, -3.066479806614716e1, 2.506628277459239];109  const b = [-5.447609879822406e1, 1.615858368580409e2, -1.556989798598866e2, 6.680131188771972e1, -1.328068155288572e1];110  const c = [-7.784894002430293e-3, -3.223964580411365e-1, -2.400758277161838, -2.549732539343734, 4.374664141464968, 2.938163982698783];111  const d = [7.784695709041462e-3, 3.224671290700398e-1, 2.445134137142996, 3.754408661907416];112  const pl = 0.02425;113  let q: number, r: number;114  if (p < pl) {115    q = Math.sqrt(-2 * Math.log(p));116    return (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1);117  }118  if (p <= 1 - pl) {119    q = p - 0.5;120    r = q * q;121    return ((((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q) / (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1);122  }123  q = Math.sqrt(-2 * Math.log(1 - p));124  return -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1);125}126127/** p-valeur bilatérale de Student. */128export function tPValue(t: number, df: number): number {129  if (!Number.isFinite(t)) return NaN;130  if (df <= 0) return NaN;131  const x = df / (df + t * t);132  return betaInc(df / 2, 0.5, x);133}134135/** Quantile t (bissection sur la CDF — suffisant). */136export function tQuantile(p: number, df: number): number {137  if (df > 200) return normalQuantile(p);138  let lo = 0;139  let hi = 50;140  const target = 2 * (1 - p); // p-valeur bilatérale correspondante141  for (let i = 0; i < 80; i++) {142    const mid = (lo + hi) / 2;143    if (tPValue(mid, df) > target) lo = mid;144    else hi = mid;145  }146  return (lo + hi) / 2;147}148149/** P(χ²_df > x). */150export function chi2Sf(x: number, df: number): number {151  if (x <= 0) return 1;152  return 1 - gammaInc(df / 2, x / 2);153}154155/** P(F_{d1,d2} > f). */156export function fSf(f: number, d1: number, d2: number): number {157  if (f <= 0) return 1;158  return betaInc(d2 / 2, d1 / 2, d2 / (d2 + d1 * f));159}160161export const stars = (p: number): string => (p < 0.001 ? "***" : p < 0.01 ? "**" : p < 0.05 ? "*" : p < 0.1 ? "·" : "");162