跳轉到

6.1 點估計#

本節目標#

  • 搞清楚「估計量」到底在估什麼：母體參數 \(\theta\) 是一個固定的數字，而估計量 \(\hat{\theta}\) 是資料的函數，因此會搖擺
  • 用偏誤（bias）與變異數（variance）兩把尺衡量一個估計量好不好
  • 理解 MSE = 變異數 + 偏誤平方，這是「不偏但變異大」與「微偏但集中」取捨的完整答案
  • 驗證兩個必背結果： \(\bar{X}\) 是不偏的；樣本變異數要**除以 \(n-1\) ** 才是對 \(\sigma^2\) 的不偏估計（呼應 1.3 與 5.4）
  • 體會「好估計量」的取捨：對稱常態用平均數、右偏資料用中位數

問題：抽到一袋樣本之後，「母體平均是多少」要怎麼講？#

第 1 章我們把資料拿來算「集中趨勢」，那一章的結論是：平均數、中位數、眾數這些數字可以描述手上的資料。但統計推論要回答的是另一個問題：

母體的参数 \(\mu\) 是一個未知但固定的數字。我們抽到 \(n\) 個樣本 \(X_1, X_2, \ldots, X_n\) ，能不能用這批資料造一個數字，來「猜」 \(\mu\) ？

這個猜測值就叫估計量（estimator） \(\hat{\theta}\) ，而它對應的參數真值 \(\theta\) 叫真值。

關鍵在於： \(\hat{\theta}\) 是樣本的函數。同一個母體，抽樣抽 100 次就得到 100 個不同的 \(\hat{\theta}\) 。所以我們不能問「它對不對」，只能問兩個統計上的問題：

  1. 中心有沒有對準？ 平均而言，它會不會系統性地偏高或偏低？——這是偏誤。
  2. 它會晃多厲害？ 反覆抽樣時估計量的散布程度有多大？——這是變異數。

這兩把尺就是本章的全部內容。

歷史：1900—1930 年代，費歇（Ronald Fisher）把「猜測」變成一門數學#

在 Fisher 之前，統計學家常把估計量叫做「最接近的近似值」，討論靠的是經驗法則。Fisher 在 1922 年的論文〈On the Mathematical Foundations of the Theory of Statistics〉裡，把問題重新定義成「決策的風險」：既然估計量會搖擺，那就用一個可以計算的函數來衡量它的表現，並且把「對準」與「晃動」拆成兩個可以分別設計的部件。

這個拆解是整個現代統計的骨架：後來的最小變異量無偏估計量（MVUE）、 Neyman–Pearson 檢定、最大似然估計，全都是這兩把尺的後代。

理論：偏誤、變異數、MSE#

定義#

對估計量 \(T = T(X_1, \ldots, X_n)\) 估計參數 \(\theta\) ：

  • 偏誤（bias）： \(\operatorname{Bias}(T) = E[T] - \theta\)

  • 變異數： \(\operatorname{Var}(T) = E\left[(T - E[T])^2\right]\)

  • 均方誤差（MSE, mean squared error）： \(\operatorname{MSE}(T) = E[(T-\theta)^2] = \operatorname{Var}(T) + \operatorname{Bias}(T)^2\)

最後那個等式是本章最重要的公式：MSE = 變異數 + 偏誤平方。它把「兩個缺點」加總成一個可以比較的數字。

注意兩個關鍵詞：

  • 不偏（unbiased） 指 \(E[T] = \theta\) ，此時偏誤為 0。不偏不等於好估計——一個估計量可以完全不偏但散布得離譜。
  • 變異數越小越好，因為它直接是「猜錯多遠」的平方的平均（ \(\operatorname{Var}(T) = E[(T-\theta)^2] - \operatorname{Bias}(T)^2\) ，不偏時兩者相等）。

必背結果一： \(\bar{X}\) 是不偏的#

\(E[\bar{X}] = E\left[\frac{1}{n}\sum_{i=1}^n X_i\right] = \frac{1}{n}\sum_{i=1}^n E[X_i] = \mu\)

這一步只用到期望值的線性性，完全不需要常態假設。變異數則是：

\(\operatorname{Var}(\bar{X}) = \frac{\sigma^2}{n}\)

也就是 5.2 節中央極限定理說的那件事：樣本越大，平均數越集中在真值附近。

必背結果二：樣本變異數為什麼要除 \(n-1\) ？#

1.3 節我們定義過樣本變異數 \(S^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X})^2\) 。現在知道那個 \(n-1\) 從哪來了：它讓 \(S^2\) 成為 \(\sigma^2\) 的不偏估計。

理由可以用自由度數一遍： \(\sum_{i=1}^n (X_i - \bar{X})^2\) 這個平方和裡， \(n\) 個偏差受到約束 \(\sum_i (X_i - \bar{X}) = 0\) ，只有 \(n-1\) 個方向能自由動（5.4 節）。而 5.4 節我們已經驗證過：

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

取期望值， \(E[\chi^2_{n-1}] = n-1\) ，所以：

\(E\left[\frac{(n-1)S^2}{\sigma^2}\right] = n-1 \;\Longrightarrow\; E[S^2] = \sigma^2\)

除以 \(n-1\) 剛好抵掉自由度少了一個的事實。 若改除以 \(n\) （有些教科書叫 \(S_n^2\) ），就會系統性低估： \(E[S_n^2] = \frac{n-1}{n}\sigma^2\) ，偏誤是 \(-\sigma^2/n\) 。

動手（一）：親眼看見 \(\bar{X}\) 的抽樣分布#

母體取 \(N(50, 12^2)\) ，重複抽樣 20000 次，每次抽 \(n\) 個，看 \(\bar{X}\) 的分布怎麼隨 \(n\) 收斂。

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
mu, sigma = 50.0, 12.0

print("=== 母體 N(50, 12^2)：重複抽樣，看 x-bar 的抽樣分布如何隨 n 收斂 ===")
print(f"{'n':>4} {'mean(x-bar)':>12} {'bias':>10} {'sd(x-bar)':>10} {'sd/sigma':>9} {'RMSE':>8}")
for n in (1, 2, 5, 10, 25, 100):
    g = r.set_seed(100 + n)
    means = [r.mean(g.rnorm(n, mean=mu, sd=sigma)) for _ in range(B)]
    bias = r.mean(means) - mu
    sd = r.sd(means)
    rmse = math.sqrt(r.mean([(v - mu) ** 2 for v in means]))
    print(f"{n:>4} {r.mean(means):>12.4f} {bias:>10.4f} {sd:>10.4f} {sd/sigma:>9.4f} {rmse:>8.4f}")
print(f"理論：sd(x-bar) = sigma/sqrt(n)；n=1 時 {sigma/math.sqrt(1):.4f}，n=25 時 {sigma/math.sqrt(25):.4f}")

輸出（實際執行）：

=== 母體 N(50, 12^2)：重複抽樣，看 x-bar 的抽樣分布如何隨 n 收斂 ===
   n  mean(x-bar)       bias  sd(x-bar)  sd/sigma     RMSE
   1      50.0557     0.0557    11.9173    0.9931  11.9171
   2      50.0574     0.0574     8.5133    0.7094   8.5133
   5      49.9913    -0.0087     5.3547    0.4462   5.3546
  10      50.0076     0.0076     3.8083    0.3174   3.8082
  25      49.9802    -0.0198     2.3924    0.1994   2.3924
 100      49.9972    -0.0028     1.2108    0.1009   1.2108
理論：sd(x-bar) = sigma/sqrt(n)；n=1 時 12.0000，n=25 時 2.4000

讀法：

  • 偏誤那一欄永遠在 0 附近跳動（最大 0.0557，只有 sigma 的 0.5%）。這就是「 \(\bar{X}\) 不偏」的實證：20000 次抽樣的平均值幾乎都在 50。
  • sd 那一欄才會隨 n 系統性下降，而且精確地遵守 \(\sigma/\sqrt{n}\) ：n=1 時 11.92（≈12），n=25 時 2.39（≈2.4），n=100 時 1.21（≈1.2）。

這兩件事要分開看：偏誤不隨 n 改善（本來就是 0），但變異會隨 n 以 \(1/\sqrt{n}\) 的速度改善。這正是 6.2 節信賴區間寬度的來源。

抽樣分布的「集中程度」#

print()
print("=== n = 1, 5, 25 的 x-bar 抽樣分布直方圖 ===")
for n in (1, 5, 25):
    g = r.set_seed(200 + n)
    means = [r.mean(g.rnorm(n, mean=mu, sd=sigma)) for _ in range(B)]
    print(f"--- n = {n} (sd = {r.sd(means):.3f}) ---")
    ascii_hist(means, nb=16)

輸出（實際執行）：

=== n = 1, 5, 25 的 x-bar 抽樣分布直方圖 ===
--- n = 1 (sd = 12.019) ---
    7.11 |# 5
   12.69 |# 31
   18.26 |## 105
   23.83 |###### 374
   29.41 |############## 866
   34.98 |############################ 1677
   40.55 |############################################# 2707
   46.13 |########################################################### 3543
   51.70 |############################################################ 3610
   57.27 |################################################## 3017
   62.85 |################################### 2111
   68.42 |################### 1156
   73.99 |######### 530
   79.57 |### 188
   85.14 |# 67
   90.71 |# 13
--- n = 5 (sd = 5.374) ---
   28.68 |# 3
   31.52 |# 9
   34.36 |# 68
   37.21 |#### 285
   40.05 |########### 777
   42.89 |########################## 1767
   45.73 |############################################# 3018
   48.57 |########################################################### 3983
   51.41 |############################################################ 4068
   54.25 |############################################## 3087
   57.09 |########################### 1826
   59.93 |########### 772
   62.77 |#### 259
   65.61 |# 65
   68.45 |# 9
   71.30 |# 4
--- n = 25 (sd = 2.399) ---
   40.78 |# 5
   41.94 |# 12
   43.10 |# 83
   44.27 |#### 230
   45.43 |########### 688
   46.59 |####################### 1469
   47.76 |###################################### 2487
   48.92 |#################################################### 3421
   50.09 |############################################################ 3916
   51.25 |################################################### 3351
   52.41 |################################### 2278
   53.58 |################### 1263
   54.74 |######### 555
   55.90 |### 178
   57.07 |# 53
   58.23 |# 11

三張圖的形狀幾乎一模一樣（都是常態），但寬度差了三倍、五倍。這是「抽樣分布集中程度」最直觀的呈現：不是形狀改變，是尺度收縮。

動手（二）：為什麼要除 \(n-1\) —— 用數字說話#

print()
print("=== 樣本變異數為什麼要除 n-1？（呼應 1.3 / 5.4）===")
g = r.set_seed(7)
s2n = []
s2m1 = []
for _ in range(B):
    s = g.rnorm(20, mean=mu, sd=sigma)
    m = r.mean(s)
    ss = sum((v - m) ** 2 for v in s)
    s2n.append(ss / 20)
    s2m1.append(ss / 19)
print(f"母體 sigma^2 = {sigma**2}")
print(f"  S2_n   (除以 n=20)  : mean = {r.mean(s2n):8.4f}  bias = {r.mean(s2n)-sigma**2:8.4f}  MSE = {r.mean([(v-sigma**2)**2 for v in s2n]):9.4f}")
print(f"  S2_(n-1)(除以 n-1=19): mean = {r.mean(s2m1):8.4f}  bias = {r.mean(s2m1)-sigma**2:8.4f}  MSE = {r.mean([(v-sigma**2)**2 for v in s2m1]):9.4f}")
print(f"  理論：MSE(S2_(n-1)) = 2*sigma^4/(n-1) = {2*sigma**4/19:.4f}，MSE(S2_n) = sigma^4*(2n-1)/n^2 = {sigma**4*39/400:.4f}")
print("  → 除以 n-1 才是對 sigma^2 的不偏估計（偏差剛好抵掉 1/(n-1) 的修正）")
print("  → 但 MSE 比較是另一回事：S2_n 雖然有偏，MSE 卻比較小！")

輸出（實際執行）：

=== 樣本變異數為什麼要除 n-1？（呼應 1.3 / 5.4）===
母體 sigma^2 = 144.0
  S2_n   (除以 n=20)  : mean = 137.2614  bias =  -6.7386  MSE = 2050.6786
  S2_(n-1)(除以 n-1=19): mean = 144.4856  bias =   0.4856  MSE = 2222.1408
  理論：MSE(S2_(n-1)) = 2*sigma^4/(n-1) = 2182.7368，MSE(S2_n) = sigma^4*(2n-1)/n^2 = 2021.7600
  → 除以 n-1 才是對 sigma^2 的不偏估計（偏差剛好抵掉 1/(n-1) 的修正）
  → 但 MSE 比較是另一回事：S2_n 雖然有偏，MSE 卻比較小！

除以 \(n\) 的版本平均只有 137.26，比真值 144 低了 6.74——正好是 \(\sigma^2/n = 144/20 = 7.2\) 的量級（抽樣誤差再疊上去）。除以 \(n-1\) 的版本平均 144.49，幾乎就是真值。

有趣的一個細節： \(S^2_{n-1}\) 的 MSE（2222）反而比 \(S^2_n\) 的 MSE（2051）稍大。這不是 \(n-1\) 比較差，而是「不偏」和「MSE 小」本來就是兩回事——這正好引出下一節的取捨。

動手（三）：好估計量的取捨#

取捨一：不偏但變異大 vs 微偏但集中#

情境 (a)：母體是對稱常態 \(N(50, 12^2)\) ，估 \(\mu\) 。比較兩個候選：

  • \(\bar{X}\) ：不偏，變異小。
  • 樣本中位數：對常態母體也幾乎不偏，但變異大約多 25%。
print()
print("=== 取捨一：不偏但變異大 vs 微偏但集中 ===")
print("(a) 對稱常態母體估 mu：x-bar（不偏）vs 中位數（也幾乎不偏，但變異較大）")
n = 25
g = r.set_seed(31)
em = []
ed = []
for _ in range(B):
    s = g.rnorm(n, mean=mu, sd=sigma)
    em.append(r.mean(s))
    ed.append(r.median(s))
print(f"  n = {n}, 母體 N({mu}, {sigma}^2)")
print(f"    x-bar  : bias = {r.mean(em)-mu:>8.4f}  sd = {r.sd(em):>7.4f}  RMSE = {math.sqrt(r.mean([(v-mu)**2 for v in em])):>7.4f}")
print(f"    median : bias = {r.mean(ed)-mu:>8.4f}  sd = {r.sd(ed):>7.4f}  RMSE = {math.sqrt(r.mean([(v-mu)**2 for v in ed])):>7.4f}")
print(f"    理論 sd(median) = 1.2533*sigma/sqrt(n) = {1.2533*sigma/math.sqrt(n):.4f}")

輸出（實際執行）：

=== 取捨一：不偏但變異大 vs 微偏但集中 ===
(a) 對稱常態母體估 mu：x-bar（不偏）vs 中位數（也幾乎不偏，但變異較大）
  n = 25, 母體 N(50.0, 12.0^2)
    x-bar  : bias =  -0.0118  sd =  2.3969  RMSE =  2.3969
    median : bias =  -0.0093  sd =  2.9856  RMSE =  2.9856
    理論 sd(median) = 1.2533*sigma/sqrt(n) = 3.0079

兩者偏誤都是 0（模擬誤差），但中位數的 sd 是 2.99 vs 平均數的 2.40，理論比值 \(1.2533\) 也對上了（3.0079 vs 2.9856）。在對稱母體上，平均數比較有效（efficient）。

情境 (b)：這個對比其實上一段已經出現過了——估 \(\sigma^2\) 時， \(S^2_n\) （除以 \(n\) ）是有偏的，但它的 MSE 卻比不偏的 \(S^2_{n-1}\) 更小。我們把它推廣到不同的 \(n\) 看：

print("(b) 微偏但更集中：估 sigma^2 用 S2_n（除以 n）還是 S2_(n-1)（除以 n-1）？")
print("    n = 20，S2_n 的偏誤是 -sigma^2/n，但 MSE 反而比較小")
for nn in (5, 10, 20, 50):
    g = r.set_seed(60 + nn)
    a = []
    b = []
    for _ in range(B):
        s = g.rnorm(nn, mean=mu, sd=sigma)
        m = r.mean(s)
        ss = sum((v - m) ** 2 for v in s)
        a.append(ss / nn)
        b.append(ss / (nn - 1))
    mse_a = r.mean([(v - sigma**2) ** 2 for v in a])
    mse_b = r.mean([(v - sigma**2) ** 2 for v in b])
    print(f"    n = {nn:>3}: S2_n   bias = {r.mean(a)-sigma**2:>8.4f}  MSE = {mse_a:>10.4f}")
    print(f"           S2_(n-1) bias = {r.mean(b)-sigma**2:>8.4f}  MSE = {mse_b:>10.4f}   (MSE 較小的是 S2_n: {mse_a < mse_b})")
print("    理論：MSE(S2_(n-1)) = 2*sigma^4/(n-1)，MSE(S2_n) = sigma^4*(2n-1)/n^2")
print(f"          n=20 時分別是 {2*sigma**4/19:.4f} 與 {sigma**4*39/400:.4f} → S2_n 的 MSE 比較小")

輸出（實際執行）：

(b) 微偏但更集中：估 sigma^2 用 S2_n（除以 n）還是 S2_(n-1)（除以 n-1）？
    n = 20，S2_n 的偏誤是 -sigma^2/n，但 MSE 反而比較小
    n =   5: S2_n   bias = -29.1448  MSE =  7410.7455
           S2_(n-1) bias =  -0.4310  MSE = 10252.2604   (MSE 較小的是 S2_n: True)
    n =  10: S2_n   bias = -14.3944  MSE =  3934.4591
           S2_(n-1) bias =   0.0062  MSE =  4601.5553   (MSE 較小的是 S2_n: True)
    n =  20: S2_n   bias =  -6.9788  MSE =  2015.2108
           S2_(n-1) bias =   0.2329  MSE =  2179.0099   (MSE 較小的是 S2_n: True)
    n =  50: S2_n   bias =  -2.8913  MSE =   816.4397
           S2_(n-1) bias =  -0.0116  MSE =   841.3994   (MSE 較小的是 S2_n: True)
    理論：MSE(S2_(n-1)) = 2*sigma^4/(n-1)，MSE(S2_n) = sigma^4*(2n-1)/n^2
          n=20 時分別是 2182.7368 與 2021.7600 → S2_n 的 MSE 比較小

這是本章最值得學會的一個判斷：

  • \(S^2_{n-1}\) ：偏誤 ≈ 0（不偏），MSE = \(2\sigma^4/(n-1)\) 。
  • \(S^2_n\) ：偏誤 = \(-\sigma^2/n\) （有偏），MSE = \(\sigma^4(2n-1)/n^2\) 。

比較兩者： \(\frac{2}{n-1} - \frac{2n-1}{n^2} = \frac{2n^2 - (2n-1)(n-1)}{n^2(n-1)} = \frac{3n-1}{n^2(n-1)} > 0\) ，所以 ** \(S^2_n\) 的 MSE 永遠比較小**（雖然差距很小）。這就是「微偏但集中」的乾淨例子。

為什麼？因為偏誤只有 \(O(1/n)\) 而變異是 \(O(1/n)\) ，**偏誤平方是 \(O(1/n^2)\) **，量級小於變異的改善量。

同一個思想的更高級版本是 James–Stein 估計量（James & Stein, 1961）：同時估計很多個參數時，把每個估計值往整體平均「收縮」，用一點點偏誤換到大幅度的變異下降。Fisher 在 1930 年代就提過類似想法，但當時沒有被推廣。

但要記住：MSE 只是很多種判準之一。在統計推論裡，「不偏」常常是刻意保留的性質，因為不偏讓我們能寫出「覆蓋率恰好是 95%」這種保證；MSE 較小但覆蓋率不對的估計量，會讓 6.2 節的信賴區間失去意義。

取捨二：平均數 vs 中位數，關鍵在於「你要估什麼」#

右偏母體（指數分布，平均 20）上，估平均和估中位數是兩件不同的事。

print()
print("=== 取捨二：右偏母體上，平均數 vs 中位數誰更有效？ ===")
g = r.set_seed(11)
xs = g.rexp(20000, rate=1.0 / 20.0)
true_med = 20.0 * math.log(2)
print(f"母體（指數，平均 20）：mean = {r.mean(xs):.3f}，median = {r.median(xs):.3f}，真中位數 = {true_med:.4f}")
for n in (10, 30):
    g = r.set_seed(300 + n)
    est_mean = []
    est_med = []
    for _ in range(B):
        s = g.rexp(n, rate=1.0 / 20.0)
        est_mean.append(r.mean(s))
        est_med.append(r.median(s))
    print(f"  n = {n}:")
    print(f"    x-bar  (目標 20.0000): bias = {r.mean(est_mean)-20.0:>8.4f}  sd = {r.sd(est_mean):>7.4f}  RMSE = {math.sqrt(r.mean([(v-20.0)**2 for v in est_mean])):>7.4f}")
    print(f"    median (目標 {true_med:.4f}): bias = {r.mean(est_med)-true_med:>8.4f}  sd = {r.sd(est_med):>7.4f}  RMSE = {math.sqrt(r.mean([(v-true_med)**2 for v in est_med])):>7.4f}")
print("  → 估『平均』用 x-bar 有效（不偏且 sd 小）；估『中位數』用樣本中位數（不偏但 O(1/n) 微偏、sd 略大）")
print("  → 真正在意的是目標參數是什麼：想估平均就別用中位數去湊")

輸出（實際執行）：

=== 取捨二：右偏母體上，平均數 vs 中位數誰更有效？ ===
母體（指數，平均 20）：mean = 19.911，median = 13.794，真中位數 = 13.8629
  n = 10:
    x-bar  (目標 20.0000): bias =  -0.0335  sd =  6.2897  RMSE =  6.2897
    median (目標 13.8629): bias =   0.9922  sd =  6.1759  RMSE =  6.2549
  n = 30:
    x-bar  (目標 20.0000): bias =  -0.0107  sd =  3.6323  RMSE =  3.6322
    median (目標 13.8629): bias =   0.3168  sd =  3.6311  RMSE =  3.6448
  → 估『平均』用 x-bar 有效（不偏且 sd 小）；估『中位數』用樣本中位數（不偏但 O(1/n) 微偏、sd 略大）
  → 真正在意的是目標參數是什麼：想估平均就別用中位數去湊

看兩件事：

  1. 樣本中位數估中位數時有 \(O(1/n)\) 的微偏（n=10 偏 0.99，n=30 偏 0.32，比例大約是 1:3，正好符合 \(1/n\) ）。這是估計中位數的固有困難（要反推分位位置），但 n 大時可以忽略。
  2. 用中位數去湊一個右偏母體的平均，會系統性地低估——這是「估錯目標」的典型災難。

結論：平均數 vs 中位數之爭不是「誰比較好」，而是「你的參數是什麼」。想估平均就用平均數（不偏、有效）；想估中位數就用中位數；資料嚴重右偏又只關心「典型值」時，中位數更穩健（不受極端值影響）。

與 R 對照#

R 沒有內建的「估計量品質」函式，但同樣的模擬可以用 replicate() 一行完成：

set.seed(42)
B <- 20000
mu <- 50; sigma <- 12

cat("=== 母體 N(50, 12^2)：x-bar 的抽樣分布 ===\n")
for (n in c(1, 5, 25, 100)) {
  means <- replicate(B, mean(rnorm(n, mu, sigma)))
  cat(sprintf("n = %3d  mean(x-bar) = %8.4f  bias = %7.4f  sd = %7.4f  理論 sd = %7.4f\n",
              n, mean(means), mean(means) - mu, sd(means), sigma / sqrt(n)))
}

cat("\n=== 樣本變異數為什麼要除 n-1 ===\n")
s <- replicate(B, { x <- rnorm(20, mu, sigma); c(sum((x - mean(x))^2) / 20, var(x)) })
cat(sprintf("S2_n     (除以 n=20)  : mean = %8.4f  bias = %8.4f\n", mean(s[1,]), mean(s[1,]) - sigma^2))
cat(sprintf("S2_(n-1)(除以 n-1=19): mean = %8.4f  bias = %8.4f\n", mean(s[2,]), mean(s[2,]) - sigma^2))

cat("\n=== x-bar vs 中位數（常態母體，n = 25）===\n")
e1 <- replicate(B, mean(rnorm(25, mu, sigma)))
e2 <- replicate(B, median(rnorm(25, mu, sigma)))
cat(sprintf("x-bar  : bias = %7.4f  sd = %6.4f\n", mean(e1) - mu, sd(e1)))
cat(sprintf("median : bias = %7.4f  sd = %6.4f\n", mean(e2) - mu, sd(e2)))

cat("\n=== 指數母體：估平均 vs 估中位數 ===\n")
true_med <- 20 * log(2)
for (n in c(10, 30)) {
  a <- replicate(B, mean(rexp(n, 1/20)))
  b <- replicate(B, median(rexp(n, 1/20)))
  cat(sprintf("n = %2d  x-bar: bias = %7.4f sd = %6.4f RMSE = %6.4f | median: bias = %7.4f sd = %6.4f RMSE = %6.4f\n",
              n, mean(a) - 20, sd(a), sqrt(mean((a - 20)^2)),
              mean(b) - true_med, sd(b), sqrt(mean((b - true_med)^2))))
}

輸出（實際執行）：

=== 母體 N(50, 12^2)：x-bar 的抽樣分布 ===
n =   1  mean(x-bar) =  49.9371  bias = -0.0629  sd = 12.0831  理論 sd = 12.0000
n =   5  mean(x-bar) =  49.9671  bias = -0.0329  sd =  5.3276  理論 sd =  5.3666
n =  25  mean(x-bar) =  50.0132  bias =  0.0132  sd =  2.4022  理論 sd =  2.4000
n = 100  mean(x-bar) =  50.0005  bias =  0.0005  sd =  1.2012  理論 sd =  1.2000

=== 樣本變異數為什麼要除 n-1 ===
S2_n     (除以 n=20)  : mean = 136.4840  bias =  -7.5160
S2_(n-1)(除以 n-1=19): mean = 143.6674  bias =  -0.3326

=== x-bar vs 中位數（常態母體，n = 25）===
x-bar  : bias = -0.0054  sd = 2.4016
median : bias = -0.0205  sd = 2.9857

=== 指數母體：估平均 vs 估中位數 ===
n = 10  x-bar: bias =  0.0543 sd = 6.3310 RMSE = 6.3311 | median: bias =  1.0457 sd = 6.1676 RMSE = 6.2554
n = 30  x-bar: bias = -0.0160 sd = 3.5958 RMSE = 3.5958 | median: bias =  0.2864 sd = 3.6016 RMSE = 3.6129

R 的結果與 Python 版一致（數字略有差異是因為亂數串不同，但所有結論相同）： \(\bar{X}\) 不偏、 \(S^2_{n-1}\) 不偏、中位數估計在右偏母體上有效率略差。

小結#

  1. 估計量是樣本的函數，會搖擺；估計量的品質用偏誤與變異數衡量，合成 MSE = 變異數 + 偏誤平方。
  2. \(\bar{X}\) 是 \(\mu\) 的不偏估計量， \(\operatorname{Var}(\bar{X}) = \sigma^2/n\) ——偏誤不隨 n 改善，但變異以 \(1/\sqrt{n}\) 改善。
  3. **樣本變異數除以 \(n-1\) ** 才是對 \(\sigma^2\) 的不偏估計，這個 \(n-1\) 就是自由度少了一個（1.3 / 5.4）。
  4. 「不偏」不是唯一目標：MSE 有時由微偏的估計量達成（本章的 \(S^2_n\) 就是例子；多參數情形則有 James–Stein 收縮）。
  5. 平均數 vs 中位數之爭的答案是「你的目標參數是什麼」：估平均用平均數，估中位數用中位數。

練習#

  1. 用模擬證明：樣本眾數（mode）不是母體眾數的不偏估計量。取一個雙峰常態混合母體，估計母體眾數位置，看重複抽樣時估計量的平均落在哪裡、偏誤多大。 （提示：先在密度上做數值搜尋找母體眾數，再用網格上的最大密度點作為樣本眾數的估計。）
  2. 驗證 \(E[S^2_n] = \frac{n-1}{n}\sigma^2\) ：對不同的 \(n = 5, 10, 20, 50\) 各模擬 5000 次，畫「 \(S^2_n\) 的平均 vs \(\frac{n-1}{n}\sigma^2\) 」的對照表。
  3. 把本節情境 (b) 推廣：對不同的 \(n = 5, 10, 20, 50\) ，算出 \(S^2_n\) 與 \(S^2_{n-1}\) 的理論 MSE（ \(2\sigma^4/(n-1)\) 與 \(\sigma^4(2n-1)/n^2\) ），與模擬值對照，並算出兩者的相對差距 \(\frac{MSE(S^2_{n-1}) - MSE(S^2_n)}{MSE(S^2_{n-1})}\) 。看 n 越大，相對差距越小。 （提示：n=5 時差距約 28%，n=50 時約 3%。）
  4. 找一個「高偏誤但低變異」的估計量：對 \(N(\mu, \sigma^2)\) 估 \(\sigma\) ，比較 \(\sqrt{S^2}\) 與「把 \(S^2\) 開根號後再除上一個小於 1 的修正係數」的版本，看哪個 MSE 小。 （提示： \(\sqrt{\cdot}\) 這一步會製造 \(O(1/n)\) 的偏誤。）
  5. 討論題：若目標參數是母體的中位數，而母體是對稱常態，則中位數與平均數都是不偏估計量。哪一個更有效？用 20000 次模擬比較兩者的 sd，並用理論說明為什麼。 （提示：對稱常態的樣本中位數的漸近效率是 \(2/\pi \approx 0.637\) 倍於平均數。）