跳轉到

5.4 卡方、t 與 F 分配的來源#

本節目標#

  • 理解三個「抽樣分配」的定義式，知道它們全部由標準常態拼出來： \(\chi^2\) = 獨立標準常態的平方和； \(t\) = 正態除以 \(\sqrt{\chi^2/\text{df}}\) ； \(F\) = 兩個 \(\chi^2\) 之比
  • 用程式從 rnorm 直接「製造」這三種隨機變數，與理論密度/機率對照
  • 模擬驗證 5.1 節的陳述： \((n-1)S^2/\sigma^2 \sim \chi^2_{n-1}\)
  • 認識 Student 的 \(t\) 的歷史（Gosset, 1908, Guinness 啤酒廠）

問題：檢定用的那三張表，是從哪裡來的？#

第 6、7 章會不斷引用三個分配：卡方（變異數的檢定與區間）、 \(t\) （平均數的檢定與區間）、 \(F\) （變異數比、變異數分析）。教科書通常直接給密度函數與查表，但它們其實都是同一個母體的推論：

母體是常態 \(N(\mu, \sigma^2)\) 時，樣本的函數 \(\bar{X}\) 、 \(S^2\) 經過適當的標準化，就會落在這三個分配上。

定義式只有三行，全部從 \(Z_1, Z_2, \ldots \overset{iid}{\sim} N(0,1)\) 出發：

$$ \chi^2_k = Z_1^2 + Z_2^2 + \cdots + Z_k^2

$$ $$ t_\nu = \frac{Z}{\sqrt{\,W/\nu\,}}, \qquad Z \sim N(0,1),\; W \sim \chi^2_\nu,\; \text{兩者獨立}

$$ $$ F_{(k,\,m)} = \frac{W_1 / k}{W_2 / m}, \qquad W_1 \sim \chi^2_k,\; W_2 \sim \chi^2_m,\; \text{兩者獨立}

$$ 幾個立刻能讀出的性質（都是模擬就能驗的）：

  • \(\chi^2_k\) ：平均 \(= k\) 、變異數 \(= 2k\) （每個 \(Z^2\) 平均 1、變異數 2）；自由度越大越對稱。
  • \(t_\nu\) ：平均 0，尾巴比常態厚（分母是隨機的，偶爾很小、把商撐大）； \(\nu \to \infty\) 時 \(\sqrt{W/\nu}\to 1\) ， \(t_\nu \to N(0,1)\) 。
  • \(F_{(k,m)}\) ：平均 \(= \frac{m}{m-2}\) （ \(m>2\) ），只取正值、右偏。

為什麼是平方和？為什麼自由度是 \(n-1\) ？#

\(S^2\) 的公式裡有一個 \(\bar{X}\) ： \(\sum(X_i - \bar{X})^2\) 的 \(n\) 個偏差受到一條約束 \(\sum(X_i - \bar{X}) = 0\) ，只剩 \(n-1\) 個自由度可動。把偏差向量標準化後，可以證明這 \(n-1\) 個「自由的方向」各自是標準常態，平方和自然就是 \(\chi^2_{n-1}\) ：

$$ \frac{(n-1)S^2}{\sigma^2} \sim \chi^2_{n-1}

$$ 而常態母體還附贈 5.1 節那條神奇的獨立性： \(\bar{X} \perp S^2\) （平均落在與偏差正交的方向上）。兩條合起來，才拼得出 \(t = \frac{\bar{X}-\mu}{S/\sqrt{n}} \sim t_{n-1}\) ——分子是常態、分母是卡方、兩者還獨立，正好是 \(t\) 的定義式。

歷史：1908 — Student 的 \(t\)#

William Gosset 是 Guinness 啤酒廠聘的牛津化學畢業生，工作是用小樣本（大麥田間試驗、每批只十幾個樣本）判斷品質。當時的統計只處理大樣本（用常態近似），他發現 \(n\) 小時估出的誤差系統性偏樂觀，於是推導了小樣本下 \(\frac{\bar{X}-\mu}{S/\sqrt{n}}\) 的精確分布。因為 Guinness 禁止員工發表業務機密，他以筆名 Student 於 1908 年發表〈The Probable Error of a Mean〉，這就是「Student's \(t\) 」名稱的由來。Fisher 後來（1920 年代）補全了嚴格證明——用的正是上面 \(t = Z/\sqrt{W/\nu}\) 的幾何。

動手（一）：自己「製造」三個分配#

全部從 rnorm 拼出來，再用 r.dchisq / r.dnorm / r.df 的理論值對照（ \(t\) 那組順便和常態比尾巴）。

import math
import rstatistics as r

def ascii_hist(x, nb=20, width=60):
    lo, hi = min(x), max(x)
    step = (hi - lo) / nb
    counts = [0] * nb
    for v in x:
        counts[min(nb - 1, int((v - lo) / step))] += 1
    mx = max(counts)
    for i, c in enumerate(counts):
        if c == 0: continue
        m = lo + (i + 0.5) * step
        bar = "#" * max(1, round(width * c / mx))
        print(f"{m:8.2f} |{bar} {c}")

B = 20000

# ---- (1) 卡方：df 個獨立標準常態的平方和，自己從 rnorm 拼 ----
print("=== 卡方：從 3 個獨立標準常態的平方和造出 chisq(df=3) ===")
g = r.set_seed(21)
chi = []
for _ in range(B):
    z1, z2, z3 = g.rnorm(3)
    chi.append(z1*z1 + z2*z2 + z3*z3)
print(f"平均 = {r.mean(chi):.3f} (理論 df = 3)，變異數 = {r.var(chi):.3f} (理論 2*df = 6)")
ascii_hist(chi, nb=18)

print()
print("=== 卡方：從 10 個獨立標準常態的平方和造出 chisq(df=10) ===")
g = r.set_seed(21)
chi10 = []
for _ in range(B):
    zs = g.rnorm(10)
    chi10.append(sum(z*z for z in zs))
print(f"平均 = {r.mean(chi10):.3f} (理論 10)，變異數 = {r.var(chi10):.3f} (理論 20)")
ascii_hist(chi10, nb=18)

# ---- (2) t：正態 / sqrt(卡方/df) ----
print()
print("=== t：Z / sqrt(W/3)，Z~N(0,1)、W~chisq(3) 獨立，造出 t(3) ===")
g = r.set_seed(22)
tv = []
for _ in range(B):
    z = g.rnorm(1)[0]
    w = sum(x*x for x in g.rnorm(3))
    tv.append(z / math.sqrt(w / 3))
print(f"平均 = {r.mean(tv):.4f} (理論 0)，sd = {r.sd(tv):.4f} (理論 sqrt(df/(df-2)) = {math.sqrt(3/1):.4f})")
print("t(3) 與標準常態的尾部對照（|t| 超過 c 的比例）：")
for c in (1.96, 3.0, 5.0):
    emp = sum(1 for v in tv if abs(v) > c) / B
    print(f"  c = {c:4.2f}:  模擬 {emp:.4f}   標準常態 {2*(1-r.pnorm(c)):.4f}   t(3) {2*(1-r.pt(c, 3)):.4f}")

# ---- (3) F：兩個卡方之比 ----
print()
print("=== F：(W1/3) / (W2/7)，W1~chisq(3)、W2~chisq(7) 獨立，造出 F(3,7) ===")
g = r.set_seed(23)
fv = []
for _ in range(B):
    w1 = sum(x*x for x in g.rnorm(3))
    w2 = sum(x*x for x in g.rnorm(7))
    fv.append((w1/3) / (w2/7))
print(f"平均 = {r.mean(fv):.4f} (理論 df2/(df2-2) = {7/5:.4f})")
print(f"超過 4.35 的比例 = {sum(1 for v in fv if v > 4.35)/B:.4f}  (理論 1 - pf(4.35,3,7) = {1-r.pf(4.35, 3, 7):.4f})")

輸出（實際執行）：

=== 卡方：從 3 個獨立標準常態的平方和造出 chisq(df=3) ===
平均 = 3.005 (理論 df = 3)，變異數 = 5.973 (理論 2*df = 6)
    0.65 |############################################################ 5423
    1.94 |############################################################ 5419
    3.24 |######################################## 3607
    4.53 |########################## 2315
    5.82 |############### 1393
    7.12 |######### 821
    8.41 |##### 497
    9.70 |### 233
   11.00 |# 126
   12.29 |# 81
   13.58 |# 35
   14.88 |# 22
   16.17 |# 7
   17.46 |# 6
   18.75 |# 6
   20.05 |# 5
   21.34 |# 2
   22.63 |# 2

=== 卡方：從 10 個獨立標準常態的平方和造出 chisq(df=10) ===
平均 = 10.040 (理論 10)，變異數 = 20.064 (理論 20)
    1.74 |#### 287
    3.82 |######################### 1636
    5.90 |##################################################### 3468
    7.99 |############################################################ 3958
   10.07 |##################################################### 3512
   12.15 |########################################## 2758
   14.23 |############################# 1886
   16.31 |################## 1163
   18.39 |######### 604
   20.47 |##### 354
   22.55 |### 211
   24.63 |# 86
   26.71 |# 42
   28.79 |# 24
   30.87 |# 2
   32.95 |# 4
   35.03 |# 1
   37.11 |# 4

=== t：Z / sqrt(W/3)，Z~N(0,1)、W~chisq(3) 獨立，造出 t(3) ===
平均 = 0.0113 (理論 0)，sd = 1.6804 (理論 sqrt(df/(df-2)) = 1.7321)
t(3) 與標準常態的尾部對照（|t| 超過 c 的比例）：
  c = 1.96:  模擬 0.1411   標準常態 0.0500   t(3) 0.1449
  c = 3.00:  模擬 0.0583   標準常態 0.0027   t(3) 0.0577
  c = 5.00:  模擬 0.0158   標準常態 0.0000   t(3) 0.0154

=== F：(W1/3) / (W2/7)，W1~chisq(3)、W2~chisq(7) 獨立，造出 F(3,7) ===
平均 = 1.4018 (理論 df2/(df2-2) = 1.4000)
超過 4.35 的比例 = 0.0505  (理論 1 - pf(4.35,3,7) = 0.0499)

三個分配都驗上了：卡方的平均/變異數對理論值、 \(t(3)\) 尾部遠厚於常態（ \(|t|>3\) 的機率是常態的 20 倍——這就是 Gosset 當年「小樣本誤差偏樂觀」的定量內容）、 \(F(3,7)\) 的平均與上 5% 臨界值都吻合。 \(t\) 的 sd 是 1.68 而理論 1.73——因為 \(t(3)\) 的變異數「存在但矩尾很重」，2 萬次模擬收斂慢，這本身也是重尾的教訓。

動手（二）：驗證 \((n-1)S^2/\sigma^2 \sim \chi^2_{n-1}\)#

把 5.1 節「留待驗證」的支票兌現：直接對常態母體重複抽樣，算統計量，看它的分布是不是 \(\chi^2_{n-1}\) 。

import math
import rstatistics as r

def ascii_hist(x, nb=20, width=60):
    lo, hi = min(x), max(x)
    step = (hi - lo) / nb
    counts = [0] * nb
    for v in x:
        counts[min(nb - 1, int((v - lo) / step))] += 1
    mx = max(counts)
    for i, c in enumerate(counts):
        if c == 0: continue
        m = lo + (i + 0.5) * step
        bar = "#" * max(1, round(width * c / mx))
        print(f"{m:8.2f} |{bar} {c}")

B = 20000

print("=== 驗證 (n-1)S^2/sigma^2 ~ chisq(n-1)：母體 N(50,12^2)，n=10 ===")
g = r.set_seed(24)
mu, sigma, n = 50.0, 12.0, 10
stat = []
for _ in range(B):
    s = g.rnorm(n, mean=mu, sd=sigma)
    stat.append((n-1) * r.var(s) / sigma**2)
print(f"平均 = {r.mean(stat):.4f} (理論 n-1 = 9)")
print(f"落在 [qchisq(0.025,9), qchisq(0.975,9)] = [{r.qchisq(0.025, 9):.4f}, {r.qchisq(0.975, 9):.4f}] 的比例 = {sum(1 for v in stat if r.qchisq(0.025,9) <= v <= r.qchisq(0.975,9))/B:.4f} (理論 0.95)")
ascii_hist(stat, nb=16)

輸出（實際執行）：

=== 驗證 (n-1)S^2/sigma^2 ~ chisq(n-1)：母體 N(50,12^2)，n=10 ===
平均 = 9.0421 (理論 n-1 = 9)
落在 [qchisq(0.025,9), qchisq(0.975,9)] = [2.7004, 19.0228] 的比例 = 0.9479 (理論 0.95)
    1.63 |####### 551
    3.87 |##################################### 2749
    6.12 |############################################################ 4412
    8.37 |########################################################### 4352
   10.62 |############################################# 3335
   12.86 |############################# 2115
   15.11 |################ 1205
   17.36 |######### 646
   19.61 |##### 333
   21.86 |## 168
   24.10 |# 74
   26.35 |# 31
   28.60 |# 17
   30.85 |# 6
   33.10 |# 4
   35.34 |# 2

平均 9.04（理論 9）、95% 區間覆蓋率 0.9479（理論 0.95）、直方圖形狀與 \(\chi^2_9\) 的右偏密度一致。第 6 章變異數的信賴區間公式

$$ \left(\frac{(n-1)S^2}{\chi^2_{\alpha/2,\,n-1}},\; \frac{(n-1)S^2}{\chi^2_{1-\alpha/2,\,n-1}}\right)

$$ 之所以是這兩個卡方分位數，就是因為這個統計量落在 \(\chi^2_{n-1}\) 上——你剛剛親手驗證了它。

與 R 對照#

R 內建 rchisq / rt / rf，可以和「自己拼」互相比對：

set.seed(22)
B <- 20000
z  <- rnorm(B);  w3 <- rchisq(B, 3);  w7 <- rchisq(B, 7)
tv <- z / sqrt(w3 / 3)
fv <- (w3 / 3) / (w7 / 7)
c(mean_tv = mean(tv), sd_tv = sd(tv))
c(mean_fv = mean(fv), upper_5 = quantile(fv, 0.95))

輸出（實際執行）：

 mean_tv        sd_tv 
-0.006854774  1.645607127 
    mean_fv upper_5.95% 
   1.395255    4.326533

（rstatistics 也提供 g.rchisq / g.rt / g.rf，只是本章刻意全部從 rnorm 拼裝，讓「來源」看得見。兩者結果在同 seed 下可比。）

小結#

  1. 三個抽樣分配一個來源： \(\chi^2\) 是常態平方和、 \(t\) 是常態除以隨機尺度、 \(F\) 是兩個卡方之比。
  2. \(\frac{(n-1)S^2}{\sigma^2} \sim \chi^2_{n-1}\) 與 \(\bar{X} \perp S^2\) 合起來推出 \(t = \frac{\bar{X}-\mu}{S/\sqrt{n}} \sim t_{n-1}\) ——第 7 章 \(t\) 檢定的引擎。
  3. \(t\) 的厚尾巴不是小毛病： \(n\) 小時用常態近似會低估誤差，Gosset（1908，Guinness）因此導出 \(t\) 分配。
  4. 「定義式 + 模擬 + 理論值對照」是本章反覆使用的三段式驗證法，之後遇到任何新分配都可以如法炮製。

練習#

  1. 從 rnorm 拼 \(t_1\) ，畫直方圖，並與 Cauchy 分配（g.rcauchy）比對——它們其實是同一個分配。模擬 \(t_1\) 的平均與 sd 會有什麼怪現象？ （提示：Cauchy 期望值不存在；模擬平均不會收斂，sd 會爆。）
  2. 驗證 \(t_\nu \to N(0,1)\) ：對 \(\nu = 1, 3, 10, 100\) 各模擬，比較 \(|T| > 1.96\) 的比例與 0.05。
  3. 母體改成非常態（如指數），重算 \((n-1)S^2/\sigma^2\) 的直方圖——還是 \(\chi^2_{n-1}\) 嗎？ （提示：不是； \(\chi^2\) 陳述依賴常態性，這正是第 7 章要談檢定穩健性的原因。）
  4. 用模擬檢查 5.1 節的另一條陳述：對每組樣本算 r.cor 所需的成對值 \((\bar{X}, S^2)\) ，畫散佈圖——常態母體下應看不出任何關聯，再換成指數母體看它出現正相關。
  5. 從 rnorm 拼 \(F(1, m)\) ，再對同一批 \(W_2\) 畫 \(\sqrt{F(1,m)}\) 的直方圖——它應該是 \(|t_m|\) 。用一小段程式確認。