橢圓的真身：拿 C 當一把尺

① 這關在解什麼問題

上一關給了你什麼 機率統計 04 ⑧ 已經把那圈橢圓畫出來了 —— 它是二維常態的等高線，也就是 03 關那張 C 的形狀被畫成地形圖。 但到機率統計 04 關結束，橢圓還只是一張圖：它告訴你哪裡人多、哪裡人少。 這一關把那張圖翻過來當工具用：同一顆橢圓，改讀成「同一圈上的每一點，離中心一樣遠」—— 橢圓就從一張圖變成一把尺。
關鍵動作只有一個：把 C 倒過來寫成 C⁻¹。 05 ④ 用 C 去找「哪個方向最散」， 這一關用 C⁻¹ 去除掉那個散開 —— 同一組主軸，反過來用。

先講一個很具體的麻煩。你在做刷卡風控，手上進來一筆交易：

你手上有這個持卡人過去兩週的紀錄，於是你做了每個工程師都會先做的事 —— 各欄各算一次 z-score （z = (值 − 平均) ÷ 標準差，02 關那套）：

• 只看金額：這個人上週在家電行刷過 5193、在電商刷過 4983。5400 的 z-score 只有 +1.45。不稀奇。
• 只看時刻：這個人是夜貓子，凌晨 02:40、03:20 都有便利商店紀錄。03:10 的 z-score 只有 −1.27。也不稀奇。

兩欄各自都放行，可是你一看就覺得毛毛的 —— 因為這個人的消費有個很清楚的規律：越晚金額越小 （深夜是買宵夜的時段）。「深夜」和「大額」單獨都正常，合起來違反規律。

這就是問題的核心：

你可能會想：那我把兩個 z-score 湊成一個向量，取長度不就好了？ 這關的 ⑦ 會用真的數字打你的臉 —— 那筆真異常算出 1.9279，而一筆完全正常的家電行消費算出 1.9277。 差 0.00013，這把尺分不出來。

這關給你的工具叫 馬氏距離（Mahalanobis distance）：一把會自動考慮「這兩欄本來就會一起動」的尺。 它的等高線不是圓，是橢圓 —— 而那顆橢圓，就是 05 關那顆虛線橢圓的真身。

② 橢圓是怎麼從 C 長出來的

這是本關最重要的一節，我會走得很慢。先把兩個字講定。

先講「等高線」

你看過地形圖。地形圖上那些一圈一圈的線，叫等高線 —— 同一條線上的每一點，海拔一樣高。 線密就是陡，線疏就是緩。你不用看整座山的立體模型，看那幾條圈就知道山長什麼樣。

我們要做一模一樣的事，只是把「海拔」換成「離資料中心多遠」。 所以下面每一張圖上的圈，意思都是：同一條圈上的每一點，在這把尺底下距離一樣遠。 圈的形狀，就是那把尺的形狀。

再把字母講定

把第三步講透：為什麼是 C⁻¹

拿 step 3 那顆紅點來手算一遍。它在 x − μ = (−1.1, 1.1)， C = [[4, 3.6], [3.6, 4]]。

先看錯的做法。兩欄的 σ 都是 √4 = 2，所以 z-score 是 (−0.55, +0.55)，取長度：

再看對的做法。先把 C 反過來（2×2 反矩陣：對調對角線、非對角線變號、整個除以 det C = 4·4 − 3.6·3.6 = 3.04）：

然後把 (x−μ)ᵀ C⁻¹ (x−μ) 這個式子整個攤開，一項一項算 —— 這就是 05 關那個二次型，只是矩陣換成 C⁻¹：

所以結論式長這樣

而這條曲線的幾何已經全部寫在 C 的特徵分解裡了 —— 這是 04、05 兩關的成果直接套用：

③ 主實驗：三個滑桿，把橢圓捏出來

上一節是我算給你看，這一節換你捏。三個滑桿直接控 C 的三個數字（C 是對稱的，所以只有三個自由度）， 橢圓即時跟著變。

④ 把這把尺命名：馬氏距離

剛剛那個式子有名字，叫馬氏距離（Mahalanobis distance，1936 年，印度統計學家 P. C. Mahalanobis）。 完整長這樣：

逐個零件拆：

寫成程式只有三行 —— 兩行取反矩陣、一行夾二次型：

// x、mu 都是 {x, y}；C 是 [[vx, cov], [cov, vy]]
function maha2(x, mu, C) {
  const Cinv = LAB.M.inv(C);          // C⁻¹：反矩陣。C 奇異時回傳 null
  const diff = LAB.M.sub(x, mu);      // x − μ：先搬到重心
  return LAB.M.quadForm(Cinv, diff);  // (x−μ)ᵀ C⁻¹ (x−μ)，回傳 d²（純量）
}
const d = Math.sqrt(maha2(x, mu, C)); // 要「幾個 σ」就開根號

比較一下歐幾里得距離的程式，你會看到唯一的差別就是中間那個 C⁻¹：

// 歐氏：d² = (x−μ)ᵀ I (x−μ)  ← 中間夾的是單位矩陣 I
function euclid2(x, mu) {
  const diff = LAB.M.sub(x, mu);
  return LAB.M.dot(diff, diff);       // 等價於 quadForm([[1,0],[0,1]], diff)
}

對照實驗：兩把尺會吵架

下面是本關的高潮。兩顆可拖的點 P 與 Q， 我同時用兩把尺量它們，然後告訴你「哪一顆比較遠」。兩把尺會給出相反的答案。

順手收一個好處：馬氏距離不怕你換單位

歐氏距離有個很煩的毛病：它的答案跟你用什麼單位有關。 把金額從「元」改記成「千元」，數字一個都沒變、資料一模一樣，歐氏距離卻整個換了一組值。

馬氏距離沒有這個問題，因為你縮放某一欄的時候，C 會跟著縮放同樣的倍數， C⁻¹ 剛好把它抵掉。⑦ 那 15 筆刷卡資料我兩種算法都跑過：

⑤ 為什麼是橢圓，而不是別的形狀

前面我都是「先給你橢圓，再解釋它」。這一節回答更根本的問題： 憑什麼等高線是橢圓？為什麼不是方形、不是菱形？

答案要從多元常態分佈（multivariate normal）借一句話。它的機率密度公式長這樣：

看清楚了嗎 —— x 只透過 d² 這一個數字影響密度，沒有別的入口。所以：

那每一圈裝了多少資料？（這裡很多人會錯）

先把「機率是什麼」升一維。機率統計 02 關講的是一維： 密度曲線底下的面積才是機率，縱軸那個高度本身不是機率。 二維的密度是一座山（機率統計 04 關那張等高線圖底下的那座山）， 所以同一句話要跟著升一維：機率 ＝ 那座山底下那一塊的體積。

於是「1σ 橢圓裝了多少機率」問的不是那條橢圓線有多長， 而是以那顆橢圓為底、垂直切下來的那一塊山有多少體積。 三維就是切一顆橢球、算超體積；再高的維度你畫不出來，但問法一字不變。

而體積會隨著維度一路漏掉。理由用一句話講得完：要落在 1σ 那一圈裡面， 每一個方向都得同時乖。一維只有一個條件要成立，二維要兩個一起成立，三維要三個 —— 每多一個維度就多乘一次折扣。一維那組 68 / 95 / 99.7 是背不過來的，實際數字是：

這三欄不是三個獨立的巧合，是同一個分布在三個自由度下的值 —— 那個分布叫卡方，⑥ 會把它拆開。

不用相信我，現場數給你看

下面這顆按鈕會真的生成 2000 個符合指定 C 的隨機點，然後一顆一顆算 d²、數落在各圈內的比例。 生成方法是 Box–Muller ＋ 著色（coloring），兩步都是自己寫的，程式碼在下面。

生成程式（兩步，都自己寫）

第一步：Box–Muller。把兩個 0~1 的均勻隨機數，變成兩個獨立的標準常態隨機數 （平均 0、變異數 1）。這是把「均勻」掰成「鐘形」的經典手法：

// 均勻 → 標準常態。取極座標：半徑用 √(−2 ln u₁)、角度用 2π·u₂，
// 這樣 (x, y) 就剛好是兩個獨立的 N(0, 1)。
function normalPair(rnd) {
  let u1 = rnd();                     // rnd() 回傳 [0, 1)
  if (u1 < 1e-12) u1 = 1e-12;         // 防 log(0) = −Infinity
  const u2 = rnd();
  const r = Math.sqrt(-2 * Math.log(u1));
  const t = 2 * Math.PI * u2;
  return { x: r * Math.cos(t), y: r * Math.sin(t) };
}

第二步：著色（coloring）。現在手上是一團標準的圓形點雲 z， 我們要把它變形成指定的 C。只要找到一個矩陣 A 滿足 A·Aᵀ = C，那麼 A·z 的共變異數矩陣就正好是 C。

找 A 的方法有兩種：Cholesky 分解，或者特徵分解取對稱平方根。 我選後者，因為 04、05 關已經把特徵分解玩熟了 —— 直接沿用同一套零件：

// C 的對稱平方根：A = v₁·√λ₁·v₁ᵀ + v₂·√λ₂·v₂ᵀ
// 直覺：沿著 v₁ 拉伸 √λ₁ 倍、沿著 v₂ 拉伸 √λ₂ 倍。圓被拉成橢圓，就這樣。
function sqrtSym(C) {
  const pr = LAB.M.eig2(C).pairs;                    // 已按 λ 由大到小
  const s1 = Math.sqrt(Math.max(0, pr[0].value));    // √λ₁
  const s2 = Math.sqrt(Math.max(0, pr[1].value));    // √λ₂
  const a = pr[0].vector, b = pr[1].vector;          // v₁, v₂（已單位化、互相垂直）
  return [
    [s1*a.x*a.x + s2*b.x*b.x,  s1*a.x*a.y + s2*b.x*b.y],
    [s1*a.x*a.y + s2*b.x*b.y,  s1*a.y*a.y + s2*b.y*b.y],
  ];
}

// 生成一顆符合 (μ, C) 的隨機點
const A = sqrtSym(C);
const z = normalPair(rnd);                           // 圓形點雲
const pt = LAB.M.add(mu, LAB.M.matVec(A, z));        // 被拉成橢圓形點雲

隨機數用的是 mulberry32（32 位元的小型 PRNG，種子固定 → 每次重新整理拿到同一批點， 數字才對得起來）。按「重新抽 2000 顆」會換種子。

⑥ d² 的分布叫卡方：門檻是從這裡查出來的

⑤ 給了你一條很好用的閉式解 P(d² ≤ k²) = 1 − exp(−k²/2)。 問題是它只有二維成立。你的資料有三欄、十欄的時候呢？ 而且下一節 ⑦ 的門檻選單會冒出 5.991、9.210 這種數字 —— 它們是從哪張表查來的？

兩個問題同一個答案：d² 這個量本身有一個分布，它叫卡方分布 （chi-squared，希臘字母 χ 唸「kai」）。它的定義短到有點失望：

為什麼 d² 剛好就是這個東西

因為 ⑤ 那個著色矩陣 A 是可以倒著走的。當時我們用 x = μ + A·z 把圓形的標準點雲拉成符合 C 的橢圓形點雲； 反過來讀就是 z = A⁻¹(x − μ) —— A⁻¹ 把橢圓揉回正圓。

而 d² 算的，正好就是「揉回正圓之後那個點離原點多遠」的平方：

門檻是把分布「反過來查」

這個動作 機率統計 02 關的累積分布與分位數那一節講過： CDF 問的是「小於等於某個值的機率是多少」，分位數是反過來問「機率 95% 的那個值是多少」。 所謂的「查卡方表」，就是在查卡方分布的分位數。

下一節 ⑦ 那個選單裡的四個數字，全部是 df = 2 的卡方分位：

⑦ 拿去用：15 筆刷卡紀錄，抓出可疑的那兩筆

回到 ① 的案子。這是某個持卡人兩週的 15 筆紀錄，兩欄是時刻與金額。 合成資料，不是真人。規律是「越晚金額越大」（下班後才買貴的、深夜只買宵夜）。 想換成你的資料，改檔案裡 CARD 那個陣列就行。

門檻要訂在哪？

如果你的資料真的是多元常態，那 d² 服從卡方分佈（chi-squared）， 自由度 = 欄位數。二維就是 2 個自由度，於是「要抓最極端的 5%」直接查表得到 d² > 5.991。 上面那個下拉選單就是這張表的四個常用值 —— ⑥ 已經把這張表拆開給你看過了， 連「為什麼 d² 會服從卡方」都有現場的曲線與模擬對帳。

二維剛好有閉式解，所以你連查表都不用：由 P(d² ≤ k²) = 1 − exp(−k²/2) 反解，門檻 = −2·ln α。 代 α = 0.05 → 5.991。代 α = 0.01 → 9.210。 這條閉式解的來歷（df = 2 的卡方剛好是指數分布）在 ⑥。

⑧ 工程師的「原來如此」

⑨ 踩坑

⑩ 這關回答了什麼