跳轉到

6.3 比例的信賴區間#

本節目標#

  • 推導最簡單的 Wald 區間 \(\hat{p} \pm z\sqrt{\hat{p}(1-\hat{p})/n}\) ，並指出它在小樣本或極端比例時的兩個致命問題
  • 學習 Wilson 區間：把「 \(p\) 的隨機性」正確地代回解聯立方程，得到永遠落在 \([0,1]\) 內的區間
  • 用 r.prop_test 產出並解讀比例的信賴區間
  • 用覆蓋率模擬量化比較 Wald 與 Wilson： \(n\) 小或 \(p\) 極端時，Wald 的覆蓋率會低到離譜
  • 學會樣本數的設計：要在 \(\pm 3\%\) 內，需要抽多少個？

問題：比例的「平均值」有沒有信賴區間？#

行銷部門想知道「我們的顧客中有多少比例會買升級方案」。抽了 \(n\) 個人，其中 \(x\) 個會買，樣本比例是：

\(\hat{p} = \frac{x}{n}\)

要估的是母體比例 \(p\) 。這看起來比平均數簡單（就是個比例），但信賴區間的問題卻更棘手，原因是： \(p\) 不是連續變數（ \(\hat{p}\) 只能取 \(0, 1/n, 2/n, \ldots, 1\) ），而且 ** \(p\) 會出現在標準誤的裡面**（因為 \(\operatorname{Var}(\hat{p}) = p(1-p)/n\) 裡的 \(p\) 是真值，不是 \(\hat{p}\) ）。

這兩點正是 Wald 區間失準的根源。

歷史：1920 — 1922 — Wald 的區間，與 Wilson 的修正#

  • 1920：Wald 發表了比例的區間估計（就是我們要推導的那個 \(`\hat{p}\pm\) 形式）。它簡單、直覺上 appealing，因此在半個世紀裡是教科書標準。
  • 1922：E. S. Wilson 在《Journal of the American Statistical Association》發表了現在稱為 Wilson score interval 的方法，處理了「 \(p\) 出現在標準誤裡」的問題。但 Wilson 當時的區間被推廣得很慢——統計軟體在往後數十年裡還是把 Wald 當預設。
  • 1998：Robert Newcombe 發表了 score 法的等價變形（後來在統計軟體裡叫 Newcombe interval），加上 1998 年 Brown–Cai–DasGupta 的 "Wilson score interval with continuity correction"，讓 Wilson 區間終於被廣泛採用。今天 R 的 prop.test 與多數軟體都給 Wilson 或其變形。

本章會用模擬把「為什麼 Wald 在小樣本會失準」完整展示出來。

理論：從抽樣分布到 Wald 區間#

抽樣分布#

\(X \sim \text{Binomial}(n, p)\) ，則：

\(E[\hat{p}] = p\) ， \(\operatorname{Var}(\hat{p}) = \frac{p(1-p)}{n}\)

由 CLT（$n` 大時）：

\(`\hat{p} \approx N\left(p, \frac{p(1-p)}{n}\right)\)

Wald 區間#

反解 \(P\left(-z_{\alpha/2} \le \frac{\hat{p}-p}{\sqrt{p(1-p)/n}} \le z_{\alpha/2}\right) = 1-\alpha\) ，把分母的 \(p\) 換成估計值 \(\hat{p}\) ：

\(P\left(\hat{p} - z_{\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \le p \le \hat{p} + z_{\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\right) \approx 1-\alpha\)

得到：

$$ \hat{p} \pm z_{\alpha/2}\,\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}

$$ 這個公式有兩個隱藏的問題：

  1. 把 \(p\) 換成 \(\hat{p}\) 這一步不對。嚴格來說，在推導中分母應該是**真值 \(p\) **，而且不等式兩邊都含 \(p\) （左邊的 \(p\) 在根號內、右邊是 \(p\) 本身），不是簡單的對稱形式。Wald 只是「換掉再對稱化」。
  2. 區間可能超出 \([0,1]\) 。若 \(\hat{p}=0\) ，標準誤是 0，區間變成 \([0,0]\) ——完全沒有資訊，但真實的 \(p\) 可能是 3%。若 \(\hat{p}=1/n\) 很小，區間下限可能是負的。

Wilson 區間：把「 \(p\) 的隨機性」放回來#

正確的做法是回到原始的聯合機率。令失誤機率 \(p\) 落在區間 \([p_L, p_U]\) 的機率為：

\(P(p_L \le p \le p_U \mid \hat{p}) = P\left(p_L \le X \le p_U \mid \hat{p} = x\right)\)

我們要求這個機率等於 \(1-\alpha\) 。在二項分布的正態近似下：

\(P(p_L \le X \le p_U) \approx \Phi\left(\frac{x + 0.5 - np_U}{\sqrt{np_U(1-p_U)}}\right) - \Phi\left(\frac{x + 0.5 - np_L}{\sqrt{np_L(1-p_L)}}\right)\)

（ \(+0.5\) 是連續性校正，先略去。）這是一個** \(p_L, p_U\) 的聯立非線性方程，沒有解析解——但可以數值求解。Wilson 的貢獻是找到一個精確到 \(O(1/n)\) ** 的閉式解：

$$ \frac{\hat{p} + z^2/(2n) \pm z\sqrt{\hat{p}(1-\hat{p})/n + z^2/(4n^2)}}{1 + z^2/n}

$$ 其中 \(z = z_{\alpha/2}\) 。

看這個公式的漂亮之處： \(\hat{p}\) 只出現在根號內部（不是作為中心），而且外面的分母 \(1+z^2/n\) 是「 \(n\) 越小、分母越靠近 0」——於是當 \(\hat{p}=0\) 時：

\(p_L = 0\) ， \(p_U = \frac{z^2/(2n) + z\sqrt{z^2/(4n^2)}}{1+z^2/n} = \frac{z^2/n}{1+z^2/n}\)

對 \(n=30, z=1.96\) 這是 \(0.1135\) —— **區間不會退化，而且永遠不超出 \([0,1]\) **。

另一個常見的替代品是 Agresti–Coull 區間（把 \(z^2/2\) 個「假成功」和「假失敗」加進分子分母），它更簡單但比 Wilson 稍寬。

動手（一）：同一組資料，兩種區間#

import math
import rstatistics as r

z975 = r.qnorm(0.975)
print(f"z_(0.975) = {z975:.6f}\n")


def wald(x, n):
    ph = x / n
    se = math.sqrt(ph * (1 - ph) / n)
    return ph - z975 * se, ph + z975 * se


def wilson(x, n):
    ph = x / n
    z2 = z975 * z975
    den = 1 + z2 / n
    cen = (ph + z2 / (2 * n)) / den
    half = z975 * math.sqrt(ph * (1 - ph) / n + z2 / (4 * n * n)) / den
    return cen - half, cen + half


print("=== 1. Wald 與 Wilson 區間公式（同一組資料）===")
cases = [
    (40, 100, 0.40),
    (5, 100, 0.05),
    (0, 30, 0.00),
    (1, 30, 0.0333),
    (30, 30, 1.00),
    (15, 30, 0.50),
    (480, 1000, 0.48),
]
print(f"{'x':>5} {'n':>6} {'p-hat':>7} {'Wald':>22} {'Wilson':>22}")
for x, n, p in cases:
    w = wald(x, n)
    wi = wilson(x, n)
    print(f"{x:>5} {n:>6} {p:>7.4f} [{w[0]:>9.4f}, {w[1]:>9.4f}] [{wi[0]:>9.4f}, {wi[1]:>9.4f}]")
print("→ Wald 在 x=0 或 x=n 時會給出退化的 [0, 0] 或 [1, 1]；Wilson 永遠落在 [0,1] 且保有寬度")

輸出（實際執行）：

z_(0.975) = 1.959964

=== 1. Wald 與 Wilson 區間公式（同一組資料）===
    x      n   p-hat                   Wald                 Wilson
   40    100  0.4000 [   0.3040,    0.4960] [   0.3094,    0.4980]
    5    100  0.0500 [   0.0073,    0.0927] [   0.0215,    0.1118]
    0     30  0.0000 [   0.0000,    0.0000] [   0.0000,    0.1135]
    1     30  0.0333 [  -0.0309,    0.0976] [   0.0059,    0.1667]
   30     30  1.0000 [   1.0000,    1.0000] [   0.8865,    1.0000]
   15     30  0.5000 [   0.3211,    0.6789] [   0.3315,    0.6685]
  480   1000  0.4800 [   0.4490,    0.5110] [   0.4492,    0.5110]
→ Wald 在 x=0 或 x=n 時會給出退化的 [0, 0] 或 [1, 1]；Wilson 永遠落在 [0,1] 且保有寬度

逐案解讀：

  • x=40, n=100（最典型的情況）：Wald [0.3040, 0.4960]，Wilson [0.3094, 0.4980]。兩者幾乎一樣，n 大時差異可以忽略。
  • x=5, n=100：Wald [0.0073, 0.0927]，Wilson [0.0215, 0.1118]。Wilson 明顯右移並加寬——因為 \(\hat{p}=0.05\) 很小， \(\hat{p}(1-\hat{p})/n = 0.000475\) 讓 Wald 嚴重低估了下限的真實不確定性。
  • x=0, n=30：Wald 給出 [0, 0]，這是完全錯誤的（我們絕對不能說「比例是 0」）；Wilson 給 [0, 0.1135]，這是合理的。
  • x=1, n=30：Wald 下限是 -0.0309（負數！）——比例不可能是負的，這是個明確的錯誤；Wilson 給 [0.0059, 0.1667]，既不越界又留了餘裕。
  • x=30, n=30：Wald 給 [1, 1]（同樣荒謬）；Wilson 給 [0.8865, 1]。
  • x=480, n=1000：兩者只差 0.0002， \(n=1000\) 時完全不用擔心。

動手（二）：r.prop_test 的區間輸出#

print()
print("=== 2. r.prop_test 的區間輸出（Wald + 連續性校正）===")
for x, n in ((40, 100), (5, 100)):
    h = r.prop_test(x, n)
    print(f"--- prop_test(x={x}, n={n}) ---")
    print(h)
    print(f"  conf_int = [{h.conf_int[0]:.7f}, {h.conf_int[1]:.7f}]，estimate = {h.estimate[0]:.7f}")
    print()

輸出（實際執行）：

=== 2. r.prop_test 的區間輸出（Wald + 連續性校正）===
--- prop_test(x=40, n=100) ---

    1-sample proportions test with continuity correction

data:  40 out of 100, null probability 0.5
X-squared = 3.61, df = 1, p-value = 0.05743
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.3047801 0.5029964
sample estimates:
  p
0.4

  conf_int = [0.3047801, 0.5029964]，estimate = 0.4000000

--- prop_test(x=5, n=100) ---

    1-sample proportions test with continuity correction

data:  5 out of 100, null probability 0.5
X-squared = 79.21, df = 1, p-value < 2.22e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.01855256 0.11829946
sample estimates:
   p
0.05

  conf_int = [0.0185526, 0.1182995]，estimate = 0.0500000

解讀 r.prop_test 的輸出：

  • data: 40 out of 100, null probability 0.5： \(x=40\) 、 \(n=100\) 、原假設的 \(p_0=0.5\) 。
  • X-squared = 3.61, df = 1：卡方統計量（配對 \(2\times 2\) 表的檢定）。
  • 95 percent confidence interval: 0.3047801 0.5029964：這是 Wald 區間加上連續性校正（R 的 prop.test 就是這樣實作的），所以它比我們手寫的純 Wald [0.3040, 0.4960] 右移了一點（0.3048 vs 0.3040，0.5030 vs 0.4960）——連續性校正把區間稍微加寬並對稱化，這正是它想達成的效果。

R 還提供另一個選擇：r.binom_test 給出 Clopper–Pearson 精確區間（直接由二項分布反解，覆蓋率保證 ≥ 95%，但偏保守、較寬）：

print()
print("=== 2b. r.binom_test 的 Clopper-Pearson 精確區間 ===")
for x, n in ((40, 100), (5, 100), (0, 30)):
    h = r.binom_test(x, n)
    print(f"--- binom_test(x={x}, n={n}) ---")
    print(h)
    print(f"  conf_int = [{h.conf_int[0]:.7f}, {h.conf_int[1]:.7f}]")
    print()

輸出（實際執行）：

=== 2b. r.binom_test 的 Clopper-Pearson 精確區間 ===
--- binom_test(x=40, n=100) ---

    Exact binomial test

data:  40 and 100
number of successes = 40, number of trials = 100, p-value = 0.05689
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.3032948 0.5027908
sample estimates:
probability of success 
                   0.4 

  conf_int = [0.3032948, 0.5027908]

--- binom_test(x=5, n=100) ---

    Exact binomial test

data:  5 and 100
number of successes = 5, number of trials = 100, p-value < 2.22e-16
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.01643188 0.11283491
sample estimates:
probability of success 
                  0.05 

  conf_int = [0.0164319, 0.1128349]

--- binom_test(x=0, n=30) ---

    Exact binomial test

data:  0 and 30
number of successes = 0, number of trials = 30, p-value = 1.863e-09
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.0000000 0.1157033
sample estimates:
probability of success 
                     0 

  conf_int = [0.0000000, 0.1157033]

三種區間在同一組資料上的對照：

x, n Wald（純） Wald + 連續性校正 Wilson Clopper–Pearson
40, 100 [0.3040, 0.4960] [0.3048, 0.5030] [0.3094, 0.4980] [0.3033, 0.5028]
5, 100 [0.0073, 0.0927] [0.0186, 0.1183] [0.0215, 0.1118] [0.0164, 0.1128]
0, 30 [0.0000, 0.0000] — [0.0000, 0.1135] [0.0000, 0.1157]

Wilson 和 Clopper–Pearson 在小樣本時幾乎一樣，這是好消息：實務上用哪一個都不會差太多。

動手（三）：覆蓋率模擬——Wald vs Wilson#

這是本章的核心實驗。對每個真實比例 \(p\) 重複抽樣 2000 次，每次算兩種區間，數一數有幾次蓋住了真值。

print()
print("=== 3. 覆蓋率比較：Wald vs Wilson ===")
B = 2000


def coverage(p_true, n, reps=B, seed=99):
    g = r.set_seed(seed)
    cw = 0
    ci = 0
    for _ in range(reps):
        x = sum(g.rbinom(n, 1, p_true))
        w = wald(x, n)
        wi = wilson(x, n)
        if w[0] <= p_true <= w[1]:
            cw += 1
        if wi[0] <= p_true <= wi[1]:
            ci += 1
    return cw / reps, ci / reps


print(f"{'p_true':>7} {'n':>6} {'Wald 覆蓋率':>12} {'Wilson 覆蓋率':>14}")
for p_true, n in ((0.5, 30), (0.5, 100), (0.5, 500),
                  (0.1, 10), (0.1, 30), (0.1, 100),
                  (0.02, 20), (0.02, 100),
                  (0.9, 20), (0.9, 30)):
    cw, ci = coverage(p_true, n)
    print(f"{p_true:>7.2f} {n:>6} {cw:>12.3f} {ci:>14.3f}")
print("（理論都應接近 0.95；Wald 在 n 小或 p 極端時會嚴重失準）")

輸出（實際執行）：

=== 3. 覆蓋率比較：Wald vs Wilson ===
 p_true      n     Wald 覆蓋率     Wilson 覆蓋率
   0.50     30        0.957          0.957
   0.50    100        0.947          0.947
   0.50    500        0.954          0.954
   0.10     10        0.653          0.933
   0.10     30        0.815          0.978
   0.10    100        0.929          0.931
   0.02     20        0.335          0.940
   0.02    100        0.869          0.948
   0.90     20        0.877          0.966
   0.90     30        0.815          0.978

這張表就是本節存在的理由：

  • \(p=0.5\) （最不極端）、 \(n \ge 30\) ：Wald 覆蓋率 0.947~0.957，完全沒問題。
  • \(p=0.1, n=10\) ：Wald 只有 0.653！意思是三分之一以上的區間根本沒蓋住真值。Wilson 有 0.933，接近 95%。
  • \(p=0.02, n=20\) ：Wald 只有 0.335，災難級。Wilson 0.940。
  • \(p=0.1, n=30\) ：Wald 0.815、Wilson 0.978（Wilson 稍微保守了一點，這是它的特徵：保證覆蓋率 ≥ 名義值）。
  • \(n=100\) 以上，兩者都回到 0.93~0.95，差異消失。

規律：

Wald 區間只有在 \(n\hat{p}\) 和 \(n(1-\hat{p})\) 都夠大（實務慣例：都 ≥ 10）時才可靠。 小樣本或極端比例時，請用 Wilson（或 Clopper–Pearson）。

動手（四）：抽樣 \(n\) 要多大？#

在做調查之前，最實際的問題是：「我要問多少人，才能把比例估在 \(\pm 3\%\) 以內？」

由 Wald 區間的半寬 \(z_{\alpha/2}\sqrt{p(1-p)/n} \le \varepsilon\) 解出：

\(n \ge \frac{z_{\alpha/2}^2\,p(1-p)}{\varepsilon^2}\)

** \(p(1-p)\) 在 \(p=0.5\) 時最大（= 0.25）**，所以在做樣本數計算時，若事先不知道 \(p\) ，一律用 0.25 最保險。

print()
print("=== 4. n 要多大才能把比例估在 ±3% 以內？ ===")
print("   （用最保守的 p-hat = 0.5，因為 p(1-p) 在 0.5 最大）")
for cl in (0.95, 0.99):
    zc = r.qnorm(1 - (1 - cl) / 2)
    for eps in (0.03, 0.05, 0.02):
        n = zc**2 * 0.25 / eps**2
        print(f"  信心水準 {cl:.2f}、誤差 ±{eps:.0%}：n >= {math.ceil(n)}")
print()
print("   若已知母體比例 p = 0.2，所需 n 更小：")
zc = z975
for eps in (0.03, 0.05):
    n = zc**2 * 0.2 * 0.8 / eps**2
    print(f"    信心水準 0.95、誤差 ±{eps:.0%}：n >= {math.ceil(n)}")

輸出（實際執行）：

=== 4. n 要多大才能把比例估在 ±3% 以內？ ===
   （用最保守的 p-hat = 0.5，因為 p(1-p) 在 0.5 最大）
  信心水準 0.95、誤差 ±3%：n >= 1068
  信心水準 0.95、誤差 ±5%：n >= 385
  信心水準 0.95、誤差 ±2%：n >= 2401
  信心水準 0.99、誤差 ±3%：n >= 1844
  信心水準 0.99、誤差 ±5%：n >= 664
  信心水準 0.99、誤差 ±2%：n >= 4147

   若已知母體比例 p = 0.2，所需 n 更小：
    信心水準 0.95、誤差 ±3%：n >= 683
    信心水準 0.95、誤差 ±5%：n >= 246

這張表的解讀：

  • ** \(\pm 3\%\) 需要 1068 人**（95% 信心水準）。這個數字常常讓人吃驚——很多人以為「1000 人調查」就夠了，其實 \(\pm 3\%\) 需要 1000 出頭。
  • 誤差不準確就要平方倍： \(\pm 5\%\) 只要 385 人（ \(3/5\) 的平方是 0.36， \(1068 \times 0.36 \approx 385\) ）， \(\pm 2\%\) 要 2401 人。誤差收窄一倍，樣本數變四倍——和 6.2 節「區間寬度減半要 \(n\) 變 4 倍」是同一條公式。
  • ** \(\pm 3\%\) 只需要 385 人（ \(\pm 3\%\) ）**……等等，讓我們看具體數字： \(\pm 5\%\) 要 385 人， \(\pm 3\%\) 要 1068 人（因為 \((5/3)^2 = 2.78\) ）， \(\pm 2\%\) 要 2401 人。
  • 若事先有把握 \(p \approx 0.2\) （ \(p(1-p) = 0.16\) 而不是 0.25），所需 n 降到 683，省了 36% 的樣本。所以「先做一個小規模預備調查」常常是划算的。

檢查：實際抽 1068 人，覆蓋率真的夠嗎？#

print()
print("=== 4b. 驗證：n = 1068 時的實際覆蓋率 ===")
g = r.set_seed(1212)
n = 1068
cw = 0
ci = 0
R = 500
for _ in range(R):
    x = sum(g.rbinom(n, 1, 0.5))
    w = wald(x, n)
    wi = wilson(x, n)
    if w[0] <= 0.5 <= w[1]:
        cw += 1
    if wi[0] <= 0.5 <= wi[1]:
        ci += 1
print(f"p_true = 0.5, n = {n}: Wald 覆蓋率 = {cw/R:.3f}，Wilson 覆蓋率 = {ci/R:.3f}")
print(f"區間半寬約 = {z975*math.sqrt(0.25/n):.4f}（目標 0.03）")

輸出（實際執行）：

=== 4b. 驗證：n = 1068 時的實際覆蓋率 ===
p_true = 0.5, n = 1068: Wald 覆蓋率 = 0.936，Wilson 覆蓋率 = 0.946
區間半寬約 = 0.0300（目標 0.03）

覆蓋率 0.936 / 0.946，都很接近 0.95。注意這裡的關鍵： \(n=1068\) 時區間的**半寬恰好是 \(\pm 3\%\) **，所以覆蓋率才會落在 95% 附近——樣本數公式把「誤差 \(\pm 3\%\) 」和「信心水準 95%」綁在一起，是同時滿足的兩件事。

動手（五）： \(n\) 變大時，兩種區間如何收斂#

print()
print("=== 5. 同一個 p-hat=0.4，n 不同時 Wald 與 Wilson 的長度 ===")
print(f"{'n':>6} {'Wald 寬度':>12} {'Wilson 寬度':>13} {'Wald 下限':>11} {'Wilson 下限':>13}")
for n in (10, 20, 50, 100, 500, 2000):
    x = round(0.4 * n)
    w = wald(x, n)
    wi = wilson(x, n)
    print(f"{n:>6} {w[1]-w[0]:>12.4f} {wi[1]-wi[0]:>13.4f} {w[0]:>11.4f} {wi[0]:>13.4f}")

輸出（實際執行）：

=== 5. 同一個 p-hat=0.4，n 不同時 Wald 與 Wilson 的長度 ===
     n      Wald 寬度     Wilson 寬度     Wald 下限     Wilson 下限
    10       0.6073        0.5191      0.0964        0.1682
    20       0.4294        0.3946      0.1853        0.2188
    50       0.2716        0.2621      0.2642        0.2761
   100       0.1920        0.1886      0.3040        0.3094
   500       0.0859        0.0856      0.3571        0.3580
  2000       0.0429        0.0429      0.3785        0.3787

** \(n \ge 100\) 之後兩者幾乎一模一樣**（0.1920 vs 0.1886、0.0429 vs 0.0429）。 \(n=10\) 時才差很多（0.6073 vs 0.5191）。

這給了一個務實的結論：

如果你的樣本數 \(n \ge 100\) （更保險的說法： \(n\hat{p} \ge 10\) 且 \(n(1-\hat{p}) \ge 10\) ），用哪個區間幾乎沒有差別；只有小樣本時才需要特別留意。

與 R 對照#

R 的 prop.test() 給出 Wald + 連續性校正 的區間（和 rstatistics 一致），binom.test() 給出 Clopper–Pearson 精確區間：

z <- qnorm(0.975)

wald <- function(x, n) {
  ph <- x/n; se <- sqrt(ph*(1-ph)/n)
  c(ph - z*se, ph + z*se)
}
wilson <- function(x, n) {
  ph <- x/n; z2 <- z^2; den <- 1 + z2/n
  cen <- (ph + z2/(2*n))/den
  half <- z*sqrt(ph*(1-ph)/n + z2/(4*n*n))/den
  c(cen - half, cen + half)
}

cat("=== Wald vs Wilson ===\n")
cases <- list(c(40,100), c(5,100), c(0,30), c(1,30), c(30,30), c(15,30), c(480,1000))
for (cs in cases) {
  x <- cs[1]; n <- cs[2]
  cat(sprintf("x=%4d n=%5d  Wald [%8.4f, %8.4f]  Wilson [%8.4f, %8.4f]\n",
              x, n, wald(x,n)[1], wald(x,n)[2], wilson(x,n)[1], wilson(x,n)[2]))
}

cat("\n=== r.prop_test 的區間輸出（Wald + 連續性校正）===\n")
print(prop.test(40, 100))
print(prop.test(5, 100))

cat("\n=== R 的 binom.test：Clopper-Pearson 精確區間 ===\n")
print(binom.test(40, 100)$conf.int)
print(binom.test(5, 100)$conf.int)
print(binom.test(0, 30)$conf.int)

cat("\n=== 覆蓋率：Wald vs Wilson（2000 次）===\n")
B <- 2000
set.seed(99)
for (pp in c(0.5, 0.1, 0.02, 0.9)) {
  for (n in c(20, 30, 100)) {
    cw <- 0; ci <- 0
    for (i in 1:B) {
      x <- sum(rbinom(n, 1, pp))
      w <- wald(x, n); wi <- wilson(x, n)
      if (w[1] <= pp && pp <= w[2]) cw <- cw + 1
      if (wi[1] <= pp && pp <= wi[2]) ci <- ci + 1
    }
    cat(sprintf("p = %.2f  n = %3d  Wald = %.3f  Wilson = %.3f\n", pp, n, cw/B, ci/B))
  }
}

cat("\n=== 要 ±3% 需要多大 n ===\n")
for (cl in c(0.95, 0.99)) {
  zc <- qnorm(1 - (1-cl)/2)
  cat(sprintf("信心水準 %.2f、±3%%：n >= %d\n", cl, ceiling(zc^2 * 0.25 / 0.03^2)))
}
cat(sprintf("已知 p = 0.2、95%%、±3%%：n >= %d\n", ceiling(z^2 * 0.2 * 0.8 / 0.03^2)))

輸出（實際執行）：

=== Wald vs Wilson ===
x=  40 n=  100  Wald [  0.3040,   0.4960]  Wilson [  0.3094,   0.4980]
x=   5 n=  100  Wald [  0.0073,   0.0927]  Wilson [  0.0215,   0.1118]
x=   0 n=   30  Wald [  0.0000,   0.0000]  Wilson [  0.0000,   0.1135]
x=   1 n=   30  Wald [ -0.0309,   0.0976]  Wilson [  0.0059,   0.1667]
x=  30 n=   30  Wald [  1.0000,   1.0000]  Wilson [  0.8865,   1.0000]
x=  15 n=   30  Wald [  0.3211,   0.6789]  Wilson [  0.3315,   0.6685]
x= 480 n= 1000  Wald [  0.4490,   0.5110]  Wilson [  0.4492,   0.5110]

=== r.prop_test 的區間輸出（Wald + 連續性校正）===

    1-sample proportions test with continuity correction

data:  40 out of 100, null probability 0.5
X-squared = 3.61, df = 1, p-value = 0.05743
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.3047801 0.5029964
sample estimates:
  p
0.4


    1-sample proportions test with continuity correction

data:  5 out of 100, null probability 0.5
X-squared = 79.21, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.01855256 0.11829946
sample estimates:
   p
0.05

=== R 的 binom.test：Clopper-Pearson 精確區間 ===
[1] 0.3032948 0.5027908
attr(,"conf.level")
[1] 0.95
[1] 0.01643188 0.11283491
attr(,"conf.level")
[1] 0.95
[1] 0.0000000 0.1157033
attr(,"conf.level")
[1] 0.95

=== 覆蓋率：Wald vs Wilson（2000 次）===
p = 0.50  n =  20  Wald = 0.957  Wilson = 0.957
p = 0.50  n =  30  Wald = 0.964  Wilson = 0.964
p = 0.50  n = 100  Wald = 0.955  Wilson = 0.955
p = 0.10  n =  20  Wald = 0.877  Wilson = 0.957
p = 0.10  n =  30  Wald = 0.824  Wilson = 0.968
p = 0.10  n = 100  Wald = 0.929  Wilson = 0.937
p = 0.02  n =  20  Wald = 0.340  Wilson = 0.935
p = 0.02  n =  30  Wald = 0.458  Wilson = 0.979
p = 0.02  n = 100  Wald = 0.869  Wilson = 0.946
p = 0.90  n =  20  Wald = 0.877  Wilson = 0.963
p = 0.90  n =  30  Wald = 0.823  Wilson = 0.975
p = 0.90  n = 100  Wald = 0.936  Wilson = 0.936

=== 要 ±3% 需要多大 n ===
信心水準 0.95、±3%：n >= 1068
信心水準 0.99、±3%：n >= 1844
已知 p = 0.2、95%、±3%：n >= 683

R 與 Python 的結果完全一致（覆蓋率數字略有差異是亂數串不同，但結論一致：Wald 在 \(p\) 極端或 \(n\) 小時覆蓋率崩潰，Wilson 穩定在 0.93~0.98）。

R 還內建了 prop.test 的 correct 參數（correct = FALSE 可關閉連續性校正）。

小結#

  1. Wald 區間 \(\hat{p} \pm z\sqrt{\hat{p}(1-\hat{p})/n}\) 來自二項分布的 CLT 近似，但把 \(p\) 換成 \(\hat{p}\) 這一步有瑕疵。
  2. Wald 的兩個毛病：(a) 區間可能超出 \([0,1]\) 或退化成一點；(b) 覆蓋率在 \(n\) 小或 \(p\) 極端時嚴重不足（ \(p=0.02, n=20\) 時只有 0.335）。
  3. Wilson 區間把 \(p\) 的隨機性正確放回聯立方程，永遠落在 \([0,1]\) 內，覆蓋率穩定在 0.93~0.98（1922 年 E. S. Wilson 提出，1998 年後才被軟體普遍採用）。
  4. r.prop_test 給 Wald + 連續性校正；r.binom_test 給 Clopper–Pearson 精確區間。
  5. 樣本數設計： \(n \ge \frac{z^2 p(1-p)}{\varepsilon^2}\) ， \(p=0.5\) 最保守。** \(\pm 3\%\) 需要 1068 人**（95% 水準）；誤差收窄一倍、 \(n\) 變四倍。

練習#

  1. 用數值求解（不用 Wilson 公式，而是直接算 \(P(p_L \le X \le p_U \mid \hat{p}=x)\) 並讓它等於 0.95）算出 \(x=5, n=100\) 的 score 區間，和 Wilson 的 [0.0215, 0.1118]$ 比較。 （提示：用 r.pbinom 或 r.dbinom 建立目標函式，然後用二分法或網格搜尋求 \(p_L, p_U\) 。）
  2. 比較 Agresti–Coull 區間與 Wilson：在覆蓋率模擬中加入 Agresti–Coull，看看它在小樣本時比 Wilson 保守多少。 （提示：Agresti–Coull 是把 \(z^2/2\) 個虛擬成功與失敗加到 \(x\) 和 \(n\) 上，再用 Wald 公式。）
  3. 找出 Wald 區間的「失效臨界點」：對 \(p=0.1\) ，讓 \(n\) 從 5 掃到 100，畫出 Wald 覆蓋率隨 \(n\) 上升的曲線，找出覆蓋率第一次達到 0.93 的 \(n\) 。 （提示：實務慣例是 \(n\hat{p} \ge 10\) ，即 \(n \ge 100\) ；但你會發現實際上需要更大。）
  4. 驗證 Clopper–Pearson 區間的覆蓋率 ≥ 95%：對 \(p=0.3\) 、 \(n=20\) 做 2000 次模擬，計算 Clopper–Pearson 的覆蓋率，看是否確實 ≥ 0.95（它應該比 Wilson 更保守）。 （提示：用 r.binom_test(x, n).conf_int 或直接用 r.qbeta 求分位數。）
  5. 討論題：市場調查中常聽到「誤差 ±3%」的宣稱。請說明為什麼這句話必須同時說明：(a) 信心水準、(b) 抽樣方式（簡單隨機 vs 分層）、(c) 是否用了保守的 \(p=0.5\) 假設。並計算若母體比例其實是 0.05， \(n=1068\) 時的實際半寬是多少。 （提示： \(p=0.05\) 時 \(p(1-p)=0.0475\) ，是 0.25 的 19%，所以實際半寬只有理論保守值的 \(\sqrt{0.19} \approx 0.436\) 倍。）