今天要解的問題

統計課本後面永遠附著幾頁密密麻麻的表：標準常態分配表、t 分配表、卡方分配表、F 分配表。學生做題目要翻到那幾頁，用手指對齊行列。

我要在網站上做這件事。兩條路：

選第二個。這篇會是全系列數學最重的一天，但實作只有大約 120 行純函式，而且完全可驗證（跟課本附表對答案）。

需要哪些函數

四個分配的 CDF，全部可以從兩個特殊函數推出來：

所以真正要寫的只有三個：erf、正規化不完全 gamma、正規化不完全 beta。

erf：Abramowitz–Stegun 有理近似

/* js/models/dist-math.js */

/* erf(x)：最大絕對誤差約 1.5e-7（A&S 7.1.26） */
function erf(x) {
  const sign = x < 0 ? -1 : 1;
  x = Math.abs(x);
  const t = 1 / (1 + 0.3275911 * x);
  const y = 1 - (((((1.061405429 * t - 1.453152027) * t) + 1.421413741) * t
                  - 0.284496736) * t + 0.254829592) * t * Math.exp(-x * x);
  return sign * y;
}

/* 標準常態 CDF */
const normCdf = z => 0.5 * (1 + erf(z / Math.SQRT2));

1.5e-7 的誤差對統計表來說綽綽有餘——課本附表只給四位小數。

（如果你需要更高精度，Cody 的 rational Chebyshev 近似可以到 1e-16，但程式碼長三倍，這裡不值得。）

不完全 gamma：卡方的核心

/* log Γ(x)：Lanczos 近似 */
function logGamma(x) {
  const g = [76.18009172947146, -86.50532032941677, 24.01409824083091,
             -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5];
  let y = x, tmp = x + 5.5;
  tmp -= (x + 0.5) * Math.log(tmp);
  let ser = 1.000000000190015;
  for (let j = 0; j < 6; j++) ser += g[j] / ++y;
  return -tmp + Math.log(2.5066282746310005 * ser / x);
}

/* 正規化不完全 gamma P(a,x)：級數展開（x 小）＋ 連分數（x 大） */
function gammaP(a, x) {
  if (x <= 0) return 0;
  if (x < a + 1) {
    /* 級數展開 */
    let ap = a, sum = 1 / a, del = sum;
    for (let n = 1; n < 300; n++) {
      ap++; del *= x / ap; sum += del;
      if (Math.abs(del) < Math.abs(sum) * 1e-15) break;
    }
    return sum * Math.exp(-x + a * Math.log(x) - logGamma(a));
  }
  /* 連分數（Lentz 演算法），算的是 Q，回傳 1-Q */
  let b = x + 1 - a, c = 1e300, d = 1 / b, h = d;
  for (let i = 1; i < 300; i++) {
    const an = -i * (i - a);
    b += 2; d = an * d + b; if (Math.abs(d) < 1e-300) d = 1e-300;
    c = b + an / c;         if (Math.abs(c) < 1e-300) c = 1e-300;
    d = 1 / d;
    const del = d * c; h *= del;
    if (Math.abs(del - 1) < 1e-15) break;
  }
  return 1 - Math.exp(-x + a * Math.log(x) - logGamma(a)) * h;
}

const chi2Cdf = (x, k) => gammaP(k / 2, x / 2);

為什麼要兩種算法：級數展開在 x < a+1 收斂快，x 大時要加幾百項還不準；連分數反之。這是數值計算的標準模式——在正確的區域用正確的展開。分界點 x < a + 1 是 Numerical Recipes 的經典建議。

logGamma 而不是 gamma：Γ(200) 會溢位成 Infinity，但 log Γ(200) 只有約 857。大數階乘一律在對數空間算，最後才 exp 回來。

不完全 beta：t 與 F 的核心

/* 正規化不完全 beta I_x(a,b) */
function betaI(a, b, x) {
  if (x <= 0) return 0;
  if (x >= 1) return 1;
  const front = Math.exp(logGamma(a + b) - logGamma(a) - logGamma(b)
                         + a * Math.log(x) + b * Math.log(1 - x));
  /* 對稱性：x > (a+1)/(a+b+2) 時換邊算，保證連分數收斂快 */
  if (x > (a + 1) / (a + b + 2)) return 1 - betaI(b, a, 1 - x);
  return front * betaCF(a, b, x) / a;
}

/* 連分數（Lentz） */
function betaCF(a, b, x) {
  const qab = a + b, qap = a + 1, qam = a - 1;
  let c = 1, d = 1 - qab * x / qap;
  if (Math.abs(d) < 1e-300) d = 1e-300;
  d = 1 / d;
  let h = d;
  for (let m = 1; m <= 300; m++) {
    const m2 = 2 * m;
    let aa = m * (b - m) * x / ((qam + m2) * (a + m2));
    d = 1 + aa * d; if (Math.abs(d) < 1e-300) d = 1e-300;
    c = 1 + aa / c; if (Math.abs(c) < 1e-300) c = 1e-300;
    d = 1 / d; h *= d * c;
    aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
    d = 1 + aa * d; if (Math.abs(d) < 1e-300) d = 1e-300;
    c = 1 + aa / c; if (Math.abs(c) < 1e-300) c = 1e-300;
    d = 1 / d;
    const del = d * c; h *= del;
    if (Math.abs(del - 1) < 1e-14) break;
  }
  return h;
}

/* t 分配 CDF（自由度 v） */
const tCdf = (t, v) => {
  const p = 0.5 * betaI(v / 2, 0.5, v / (v + t * t));
  return t > 0 ? 1 - p : p;
};

/* F 分配 CDF（自由度 d1, d2） */
const fCdf = (f, d1, d2) =>
  f <= 0 ? 0 : betaI(d1 / 2, d2 / 2, d1 * f / (d1 * f + d2));

betaI 裡那行對稱性交換是關鍵：I_x(a,b) = 1 - I_{1-x}(b,a)。連分數在 x 靠近 1 時收斂很慢，換邊之後永遠在快的那一側。少了這行，t=0.1, v=100 這種輸入會算幾百次迭代還不準。

反函數：雙向查詢

課本用表是兩個方向的：「z=1.96 對應多少機率」和「95% 對應多少 z」。第二個方向需要反函數。

不必為每個分配寫解析反函數——二分搜尋就夠：

/* 通用分位數：對單調 CDF 做二分搜尋 */
function quantile(cdf, p, lo, hi, iter = 200) {
  if (p <= 0) return lo;
  if (p >= 1) return hi;
  for (let i = 0; i < iter; i++) {
    const mid = (lo + hi) / 2;
    if (cdf(mid) < p) lo = mid; else hi = mid;
  }
  return (lo + hi) / 2;
}

const zCrit    = p          => quantile(normCdf, p, -40, 40);
const tCrit    = (p, v)     => quantile(t => tCdf(t, v), p, -1e3, 1e3);
const chi2Crit = (p, k)     => quantile(x => chi2Cdf(x, k), p, 0, 1e6);
const fCrit    = (p, d1, d2)=> quantile(f => fCdf(f, d1, d2), p, 0, 1e6);

200 次二分把區間縮小 2⁻²⁰⁰ 倍，遠超浮點精度——收斂性不用擔心，因為 CDF 單調。每次呼叫要算 200 次 CDF，聽起來很多，但實測整個查詢在 1 毫秒內。不要為了優雅去手寫 Newton 迭代（還得提供 PDF、還得處理不收斂），二分法笨但絕不失敗。

UI：雙向輸入 + 分配曲線

const TableView = {
  render(root, state) {
    const { dist, mode, x, p, df1, df2 } = state;
    root.innerHTML = `
      <div class="tbl-tabs">
        ${["z", "t", "chi2", "F"].map(k =>
          `<button class="tbl-tab ${k === dist ? "active" : ""}" data-dist="${k}">${LABEL[k]}</button>`).join("")}
      </div>
      <div class="tbl-body">
        ${dfInputs(dist, df1, df2)}
        <div class="tbl-io">
          <label>統計量 x <input id="in-x" type="number" step="0.0001" value="${x.toFixed(4)}"></label>
          <span class="tbl-arrow">⇄</span>
          <label>累積機率 P(X ≤ x) <input id="in-p" type="number" step="0.0001" min="0" max="1" value="${p.toFixed(6)}"></label>
        </div>
        <div class="tbl-derived">
          右尾 P(X &gt; x) = <b>${(1 - p).toFixed(6)}</b>
          ${dist === "z" || dist === "t" ? `· 雙尾 = <b>${(2 * Math.min(p, 1 - p)).toFixed(6)}</b>` : ""}
        </div>
        ${this.curveSVG(dist, state)}
      </div>`;
  },

⇄ 是這個頁面的核心互動：改左邊算右邊，改右邊算左邊。這正是課本查表的兩種用法，而紙本做不到即時切換。

分配曲線用 SVG 畫，把 P(X ≤ x) 的區域塗色：

  curveSVG(dist, { x, df1, df2 }) {
    const pdf = PDF[dist];                       // 對應的機率密度函數
    const [lo, hi] = RANGE[dist](df1, df2);
    const N = 240, W = 640, H = 200;
    const pts = [];
    let maxY = 0;
    for (let i = 0; i <= N; i++) {
      const t = lo + (hi - lo) * i / N;
      const y = pdf(t, df1, df2);
      maxY = Math.max(maxY, y);
      pts.push([t, y]);
    }
    const sx = t => ((t - lo) / (hi - lo)) * W;
    const sy = y => H - (y / maxY) * (H - 20);

    const line = pts.map(([t, y], i) => `${i ? "L" : "M"} ${sx(t).toFixed(1)} ${sy(y).toFixed(1)}`).join(" ");
    /* 左尾填色：曲線到 x 為止，再沿底邊收回 */
    const fillPts = pts.filter(([t]) => t <= x);
    const fill = fillPts.length > 1
      ? `M ${sx(lo)} ${H} ` + fillPts.map(([t, y]) => `L ${sx(t).toFixed(1)} ${sy(y).toFixed(1)}`).join(" ")
        + ` L ${sx(Math.min(x, hi)).toFixed(1)} ${H} Z`
      : "";

    return `<svg class="dist-svg" viewBox="0 0 ${W} ${H}" role="img"
                 aria-label="${LABEL[dist]} 分配曲線，已標示 x=${x.toFixed(4)} 左側面積">
      ${fill ? `<path class="dist-fill" d="${fill}"/>` : ""}
      <path class="dist-line" d="${line}"/>
      <line class="dist-mark" x1="${sx(x)}" y1="0" x2="${sx(x)}" y2="${H}"/>
    </svg>`;
  },

視覺化在這裡不是裝飾——「機率就是曲線下的面積」是初學者最難建立的直覺，一張會即時變動的塗色圖比十行文字有效。這也直接呼應 Day 9 的觀點：互動比文字更能傳達直覺。

驗證：跟課本附表對答案

這是今天最重要的部分。數值程式必須有基準測試，否則你根本不知道它對不對。

/* scripts/check-dist.js */
const D = require("../js/models/dist-math-node.js");   // 同一份程式，加個 module.exports 包裝

const CASES = [
  /* [說明, 算出來的值, 課本附表值, 容許誤差] */
  ["z: P(Z≤1.96)",        D.normCdf(1.96),          0.9750, 5e-5],
  ["z: P(Z≤1.645)",       D.normCdf(1.645),         0.9500, 5e-5],
  ["z: P(Z≤-2.58)",       D.normCdf(-2.58),         0.0049, 5e-5],
  ["z 反查: 97.5%",        D.zCrit(0.975),           1.9600, 5e-4],
  ["t: t(0.975, 10)",     D.tCrit(0.975, 10),       2.2281, 5e-4],
  ["t: t(0.95, 30)",      D.tCrit(0.95, 30),        1.6973, 5e-4],
  ["t: v=∞ 應趨近 z",      D.tCrit(0.975, 100000),   1.9600, 1e-3],
  ["chi2: χ²(0.95, 5)",   D.chi2Crit(0.95, 5),     11.0705, 1e-3],
  ["chi2: χ²(0.05, 20)",  D.chi2Crit(0.05, 20),    10.8508, 1e-3],
  ["F: F(0.95, 3, 12)",   D.fCrit(0.95, 3, 12),     3.4903, 1e-3],
  ["F: F(0.95, 1, 10) = t²", D.fCrit(0.95, 1, 10),
                              Math.pow(D.tCrit(0.975, 10), 2), 1e-3],
];

let bad = 0;
for (const [name, got, want, tol] of CASES) {
  const ok = Math.abs(got - want) <= tol;
  if (!ok) bad++;
  console.log(`${ok ? "✓" : "✗"} ${name}: got ${got.toFixed(6)} want ${want.toFixed(6)}`);
}
process.exit(bad ? 1 : 0);

最後兩個案例是數學恆等式，不是查表值：

• t 分配自由度趨近無限時應收斂到常態（t(0.975, ∞) = z(0.975) = 1.96）。
• F(1-α, 1, ν) = [t(1-α/2, ν)]²。

這種交叉驗證比對表更有價值：它同時驗了兩個獨立實作的一致性。如果我的 betaI 有 bug，這兩條幾乎一定會失敗。

踩到的雷

遞迴的 betaI 忘記處理邊界。 betaI(a, b, 1) 會呼叫 betaI(b, a, 0)，如果沒有 x <= 0 return 0 的守衛，就是無窮遞迴 + log(0) = -Infinity。邊界條件要在函式最前面就處理掉。

quantile 的初始上界不能亂設。 我一開始 chi2Crit 用 hi = 1000，結果自由度 500 的卡方分位數超過 1000，二分永遠貼在上界，回傳 1000。改成 1e6 就好。二分搜尋的初始區間必須真的包住答案，這是最容易忽略的失敗模式——而且它不會報錯，只會回傳邊界值。

浮點下溢。 Lentz 連分數演算法裡那些 if (Math.abs(d) < 1e-300) d = 1e-300 看起來很魔法，但沒有它們，某些參數組合會除以 0 得到 NaN。這是原始論文就有的保護，照抄別自己簡化。

同一份數學程式要能在瀏覽器與 Node 都跑。 我的做法是主檔給瀏覽器（全域函式），另外一個 20 行的包裝檔給 Node（require 主檔內容 + module.exports）。或者用 Day 3 的 vm 技巧把它讀進測試。不要維護兩份實作。

驗證

node scripts/check-dist.js       # 11 個案例全過才 exit 0
node scripts/verify.js
python3 -m http.server 8901
# http://localhost:8901/tables.html

手動：切換四個分配、輸入 x 看 p 變、輸入 p 看 x 變、改自由度看曲線變形、確認 F 分配的兩個自由度欄位都有作用。

順手接進 CI（Day 22）——數值程式一旦「重構」很容易靜默退化，這 11 個案例就是防線。

小結與明天預告

今天的重點：

1. 能算就不要存表：資料量為零、任意輸入、雙向查詢。
2. 數值計算的兩個通則：在正確區域用正確的展開（級數 vs 連分數）、大數在對數空間算。
3. 反函數用二分搜尋：笨但絕不失敗，不需要 Newton 也不需要解析式。
4. 數值程式必須有基準測試，而且要包含數學恆等式的交叉驗證。

明天做章末總測驗——從 Day 8 的「一課一題 check-in」升級成多題計分、選項亂序、錯題回顧。順便處理 Day 8 留下的坑：亂序之後 data-answer 這種位置索引就不能用了。