母比率の信頼区間（正規分布へ近似を利用するもの）

これは何か

• 母比率の信頼区間について調べたことをまとめたもの
• 信頼区間構築のやり方はいくつもあるが, 本記事では正規分布への近似を利用するものについてまとめた
  • 具体的には, 以下の3つの信頼区間について整理した
    • Wald interval
    • Wilson interval
    • Agresti–Coull interval

記号の定義

• Binom$(n, p)$: 2項分布
• $z_{\alpha}$: 標準正規分布の上側$100 \alpha$%点
• 本記事では表記簡略化のため、信頼区間は $CI = c \pm d$ の形で表記
  • $CI = [c - d, c + d]$を意味している

各信頼区間の概要

以後, $X \sim$ Binom$(n, p)$とし, $\hat{p} = \frac{X}{n}$とする.

Wald interval

• 2項分布の正規分布近似から導かれる信頼区間
• 入門的な統計学のテキストで頻繁に見かける

$$CI_{\rm{Wald}} = \hat{p} \pm z_{\alpha/2} \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}.$$

Wald intervalの導出

まず, 中心極限定理より, $n \rightarrow \infty$で$\frac{\hat{p} - p}{\sqrt{\frac{p(1-p)}{n}}}$は標準正規分布に分布収束することがわかる.

信頼区間は, 帰無仮説$H_{0}: p=p_{0}$が棄却されない$p_{0}$の範囲である. そこで上記統計量を利用し, 信頼水準$100(1-\alpha)$%の両側検定を考える. この検定で帰無仮説$H_{0}$が棄却されない範囲が信頼区間となるので, 以下の条件を満たす$p_{0}$の範囲が求める信頼区間となる:

$$\frac{\left| \hat{p} - p_{0} \right|}{\sqrt{\frac{p_{0}(1-p_{0})}{n}}} \le z_{\alpha/2}.$$

この不等式を$p_{0}$について解けば, 母比率の信頼区間が求まる. この不等式は2次不等式に帰着され, そのまま解くことはできる（そうして得られるのが Wilson interval）. ここで「分母に含まれる母比率パラメータ$p_{0}$を標本比率$\hat{p}$と置き換える」ことで, この不等式を1次不等式にして解くことで得られる信頼区間がWald intervalである.

この考え方自体は理論的に正当化できるが, 本題から外れるため, 補足情報として触れるにとどめておく.

（補足）

示すべきことは, 分母のパラメータを置き換えた統計量$\frac{\hat{p} - p_{0}}{\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}}$でも標準正規分布に分布収束することである.

まず, 注目している統計量を$\frac{\hat{p} - p_{0}}{\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}} = \frac{\hat{p} - p_{0}}{\sqrt{\frac{p_{0}(1-p_{0})}{n}}} \sqrt{\frac{p_{0}(1-p_{0})}{\hat{p}(1-\hat{p})}}$のように2つの統計量の積として表現する. 1つ目の統計量は先述の通り, 標準正規分布に分布収束する（中心極限定理）. ここで, 大数の弱法則から$\hat{p}$は$p_{0}$に確率収束することを踏まえると, 2つ目の統計量は$\sqrt{\frac{p_{0}(1-p_{0})}{p_{0}(1-p_{0})}} = 1$に確率収束することがわかる. よってSlutskyの定理より, 2つの統計量の積の分布は標準正規分布に分布収束するとわかる.

ここで本題に戻る. 「分母に含まれる母比率パラメータ$p_{0}$を標本比率$\hat{p}$と置き換える」ことで, 元の不等式は次のように$p_{0}$に関する1次不等式になる:

$$\frac{\left| \hat{p} - p_{0} \right|}{\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}} \le z_{\alpha/2} \\ \Leftrightarrow - z_{\alpha/2} \le \frac{\hat{p} - p_{0}}{\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}} \le z_{\alpha/2}.$$

これを$p_{0}$について解けば

$$\hat{p} - z_{\alpha/2} \sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \le p_{0} \le \hat{p} + z_{\alpha/2} \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$$

となり, Wald interval

$$CI_{\rm{Wald}} = \hat{p} \pm z_{\alpha/2} \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$$

が得られる.

Wilson interval

• Wald intervalの導出の過程で出てくる不等式をそのまま解くことで得られる信頼区間
• 一見すると複雑な式だが, 区間の中央の値は $\hat{p}$ と $\frac{1}{2}$ の加重平均となっている
  • $\frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n}}{1 + \frac{z_{\alpha/2}^{2}}{n}} = \frac{n}{n + z_{\alpha/2}^{2}} \hat{p} + \frac{z_{\alpha/2}^{2}}{n + z_{\alpha/2}^{2}} \frac{1}{2}$
  • Wald intervalの中央の値は$\hat{p}$なので, Wilson intervalの中央の値はそれと比べて$p = \frac{1}{2}$に寄った形となる

$$CI_{\rm{Wilson}} = \frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \pm \frac{z_{\alpha/2}}{\sqrt{n}} \sqrt{\hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n}}}{1 + \frac{z_{\alpha/2}^{2}}{n}}.$$

Wilson intervalの導出

基本的な考え方はWald intervalと同じだが, 途中で出てきた以下の不等式をそのまま解くのが異なる点である:

$$\frac{\left| \hat{p} - p_{0} \right|}{\sqrt{\frac{p_{0}(1-p_{0})}{n}}} \le z_{\alpha/2}.$$

両辺を2乗して分母を払うと

$$\left( \hat{p} - p_{0} \right)^{2} \le z_{\alpha/2}^{2} \frac{p_{0}(1-p_{0})}{n}$$

となり, $p_{0}$について整理すると

$$\left( 1 + \frac{z_{\alpha/2}^{2}}{n} \right) p_{0}^{2} - \left( 2 \hat{p} + \frac{z_{\alpha/2}^{2}}{n} \right) p_{0} + \hat{p}^{2} \le 0$$

という$p_{0}$についての2次不等式が得られる. あとはこれを解けばよい.

等号が成立する$p_{0}$は, 左辺=0の2次方程式を解けば得られる. 等号が成立する$p_{0}$の値は, 2次方程式の解の公式より

$$\begin{aligned} p_{0} &= \frac{2 \hat{p} + \frac{z_{\alpha/2}^{2}}{n} \pm \sqrt{\left( 2 \hat{p} + \frac{z_{\alpha/2}^{2}}{n} \right)^{2} - 4 \left( 1 + \frac{z_{\alpha/2}^{2}}{n} \right) \hat{p}^{2}}}{2 \left( 1 + \frac{z_{\alpha/2}^{2}}{n} \right)} \\ &= \frac{2 \hat{p} + \frac{z_{\alpha/2}^{2}}{n} \pm \frac{2 z_{\alpha/2}}{\sqrt{n}} \sqrt{\hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n}}}{2 \left( 1 + \frac{z_{\alpha/2}^{2}}{n} \right)} \\ &= \frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \pm \frac{z_{\alpha/2}}{\sqrt{n}} \sqrt{\hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n}}}{1 + \frac{z_{\alpha/2}^{2}}{n}} \end{aligned}$$

となる.

以上よりWilson interval

$$CI_{\rm{Wilson}} = \frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \pm \frac{z_{\alpha/2}}{\sqrt{n}} \sqrt{\hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n}}}{1 + \frac{z_{\alpha/2}^{2}}{n}}$$

が得られる.

Agresti–Coull interval

• Wald intervalをベースに$n, X, \hat{p}$を補正して得られる信頼区間
  • 区間の中央の値がWilson intervalの中央の値と同じになるように補正

$$\begin{aligned} CI_{\rm{Agresti–Coull}} &= \tilde{p} \pm z_{\alpha/2} \sqrt{\frac{\tilde{p}(1-\tilde{p})}{\tilde{n}}} \\ &= \frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \pm \frac{z_{\alpha/2}}{\sqrt{n + z_{\alpha/2}^{2}}} \sqrt{\left( \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \right) \left( 1 - \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \right)}}{1 + \frac{z_{\alpha/2}^{2}}{n}}. \end{aligned}$$

ただし, $\tilde{n}, \tilde{X}, \tilde{p}$は以下で定義される:

• $\tilde{n} = n + z_{\alpha/2}^{2}$
• $\tilde{X} = X + \frac{z_{\alpha/2}^{2}}{2}$
• $\tilde{p} = \frac{\tilde{X}}{\tilde{n}}$.

Agresti–Coull intervalの導出

Wald intervalは

$$CI_{\rm{Wald}} = \hat{p} \pm z_{\alpha/2} \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$$

であった. 式の形自体はシンプルではあるので, この式の形を維持しつつ, $n, X, \hat{p}$を補正することでより良い信頼区間を得る, ということを考える.

1つの考え方として, 中央の値がWilson intervalの中央の値 $\frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n}}{1 + \frac{z_{\alpha/2}^{2}}{n}} = \frac{X + \frac{z_{\alpha/2}^{2}}{2}}{n + z_{\alpha/2}^{2}}$になるように調整する, ということを考える.

ここで, やや天下り的ではあるが

• $\tilde{n} = n + z_{\alpha/2}^{2}$
• $\tilde{X} = X + \frac{z_{\alpha/2}^{2}}{2}$
• $\tilde{p} = \frac{\tilde{X}}{\tilde{n}}$

のように定義してみる. すると, $\tilde{p}$ がちょうどWilson intervalの中央の値に一致する.

この$\tilde{n}, \tilde{X}, \tilde{p}$を使ってWald intervalを構成したものがAgresti–Coull intervalとなる:

$$\begin{aligned} CI_{\rm{Agresti–Coull}} &= \tilde{p} \pm z_{\alpha/2} \sqrt{\frac{\tilde{p}(1-\tilde{p})}{\tilde{n}}} \\ &= \frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n}}{1 + \frac{z_{\alpha/2}^{2}}{n}} \pm z_{\alpha/2} \sqrt{\frac{1}{n + z_{\alpha/2}^{2}} \frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n}}{1 + \frac{z_{\alpha/2}^{2}}{n}} \left( 1 - \frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n}}{1 + \frac{z_{\alpha/2}^{2}}{n}} \right)} \\ &= \frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \pm \frac{z_{\alpha/2}}{\sqrt{n + z_{\alpha/2}^{2}}} \sqrt{\left( \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \right) \left( 1 - \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \right)}}{1 + \frac{z_{\alpha/2}^{2}}{n}}. \end{aligned}$$

信頼区間の観察

信頼区間の数値例

実際に$\alpha = 0.05, n = 20$のケースを例に, 3つの信頼区間の上限・下限の値をプロットしてみる.

（可視化コード）

confidence_intervals.py

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm


def wald_confidence_interval(n: int, x: int, alpha: float = 0.05) -> tuple[float, float]:
    z = norm.ppf(1 - alpha / 2)
    hat_p = x / n
    deviation = z * np.sqrt(hat_p * (1 - hat_p) / n)
    lower = hat_p - deviation
    upper = hat_p + deviation
    return (lower, upper)


def wilson_confidence_interval(n: int, x: int, alpha: float = 0.05) -> tuple[float, float]:
    z = norm.ppf(1 - alpha / 2)
    hat_p = x / n
    center = hat_p + (z**2) / (2*n)
    deviation = z / np.sqrt(n) * np.sqrt(hat_p * (1 - hat_p) + (z**2) / (4*n))
    denominator = 1 + (z**2)/n
    lower = (center - deviation) / denominator
    upper = (center + deviation) / denominator
    return (lower, upper) 


def agresti_coull_confidence_interval(n: int, x: int, alpha: float = 0.05) -> tuple[float, float]:
    z = norm.ppf(1 - alpha / 2)
    tilde_n = n + z**2
    tilde_x = x + (z**2)/2
    tilde_p = tilde_x / tilde_n
    deviation = z * np.sqrt(tilde_p * (1 - tilde_p) / tilde_n)
    lower = tilde_p - deviation
    upper = tilde_p + deviation
    return (lower, upper)


def main():
    n = 20
    alpha = 0.05

    xs = np.arange(n + 1)
    wald_lower, wald_upper = np.zeros(n + 1), np.zeros(n + 1)
    wilson_lower, wilson_upper = np.zeros(n + 1), np.zeros(n + 1)
    agresti_coull_lower, agresti_coull_upper = np.zeros(n + 1), np.zeros(n + 1)

    for x in xs:
        wald_ci = wald_confidence_interval(n, x, alpha)
        wald_lower[x], wald_upper[x] = wald_ci[0], wald_ci[1]

        wilson_ci = wilson_confidence_interval(n, x, alpha)
        wilson_lower[x], wilson_upper[x] = wilson_ci[0], wilson_ci[1]

        agresti_coull_ci = agresti_coull_confidence_interval(n, x, alpha)
        agresti_coull_lower[x], agresti_coull_upper[x] = agresti_coull_ci[0], agresti_coull_ci[1]

    plt.plot(xs, np.zeros(len(xs)), linestyle="-", color="black")
    plt.plot(xs, np.ones(len(xs)), linestyle="-", color="black")

    plt.plot(xs, wald_lower, label="wald", linestyle=":", marker="o", color="blue", markerfacecolor='none')
    plt.plot(xs, wald_upper, linestyle=":", marker="o", color="blue", markerfacecolor='none')

    plt.plot(xs, wilson_lower, label="wilson", linestyle=":", marker="x", color="red")
    plt.plot(xs, wilson_upper, linestyle=":", marker="x", color="red")

    plt.plot(xs, agresti_coull_lower, label="agresti-coull", linestyle=":", marker="+", color="green")
    plt.plot(xs, agresti_coull_upper, linestyle=":", marker="+", color="green")

    plt.xlabel("$x$")
    plt.ylabel("$p$")
    plt.title(f"Confidence Intervals (n={n} / confidence level={1-alpha})")
    plt.legend(loc="upper left")
    plt.show()


if __name__ == "__main__":
    main()

横軸を$X$の値, 縦軸を母比率の信頼区間の上限・下限の値としてプロットしている:

すると, グラフから以下のことが読み取れる:

• Wald interval は $X = 0, n$ において信頼区間が1点のみになる
• Wilson interval は0以上1以下の範囲に収まっている
  • 逆に, Wald interval・Agresti–Coull interval は0未満・1より大きい範囲を含むことがある
• Agresti–Coull interval は Wilson interval を含む

これらの性質が成り立つことは比較的簡単な計算で確認ができるため, 後ほど確認していく.

信頼区間の被覆確率

1つ前のセクションでは, 信頼区間の実現値の挙動について確認をした. その際, 「真の母比率」が信頼区間の実現値に含まれるか, といった議論はしていなかった. そこで本セクションでは, （実世界が2項分布で表現できる前提のもとでの）真の母比率と信頼区間の関係性について確認する.

被覆確率を「2項分布モデルのもとで, （確率変数としての）信頼区間が真の母比率を含む確率」と定義し, その値がどうなるかを調べてみる. 被覆確率を数式で表現すると以下のようになる（被覆確率は, 2項分布のパラメータ$n, p$と信頼区間$CI$ごとに定まる数値）:

$${\rm P}_{p} \left( CI \ni p \right) = \sum_{x: CI \ni p} \binom{x}{n} p^{x} \left( 1 - p \right)^{n-x}.$$

ここで, 3点補足を付け加えておく:

• 被覆確率の高さだけで信頼区間の良し悪しを論じるのはNG
  • 極端な話, 常に信頼区間を$\left[ 0, 1 \right]$としておけば被覆確率は1（100%）になるが, 「母比率の信頼区間は$\left[ 0, 1 \right]$です」と言われても何の情報も得られない
  • 信頼区間の幅など, 他の観点も併せて評価する必要がある
• ある程度妥当性を感じられる被覆確率の値は, $1-\alpha$（$100(1-\alpha)$%）以上
  • Wald interval の議論を思い出すと, 本記事で扱った信頼区間はいずれも「信頼水準$100(1-\alpha)$%の両側検定から信頼区間を構築する」というアイデアに基づいている
  • 仮説検定の観点から見ると, 第1種の過誤を$\alpha$（$100 \alpha$%）以下に抑えることを期待している
  • これを被覆確率の観点で言い換えれば, 被覆確率が$1-\alpha$（$100(1-\alpha)$%）以上になることを期待している, ということになる
• 確率変数としての信頼区間と, 実際のデータから計算した実現値としての信頼区間の区別は必要
  • 被覆確率の議論では, 確率変数としての信頼区間を念頭においている
  • 実際のデータから計算した実現値としての信頼区間の場合, その区間が真の母比率を含むか否かはデータを取った時点で確定しているため, 「その区間が真の母比率を含む確率」なるものを考えたところで0 or 1にしかならない

実際に$\alpha = 0.05, n = 20$のケースを例に, 3つの信頼区間の被覆確率をプロットしてみる.

（可視化コード）

coverage_probabilities.py

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm, binom


def wald_confidence_interval(n: int, x: int, alpha: float = 0.05) -> tuple[float, float]:
    z = norm.ppf(1 - alpha / 2)
    hat_p = x / n
    deviation = z * np.sqrt(hat_p * (1 - hat_p) / n)
    lower = hat_p - deviation
    upper = hat_p + deviation
    return (lower, upper)


def wilson_confidence_interval(n: int, x: int, alpha: float = 0.05) -> tuple[float, float]:
    z = norm.ppf(1 - alpha / 2)
    hat_p = x / n
    center = hat_p + (z**2) / (2*n)
    deviation = z / np.sqrt(n) * np.sqrt(hat_p * (1 - hat_p) + (z**2) / (4*n))
    denominator = 1 + (z**2)/n
    lower = (center - deviation) / denominator
    upper = (center + deviation) / denominator
    return (lower, upper) 


def agresti_coull_confidence_interval(n: int, x: int, alpha: float = 0.05) -> tuple[float, float]:
    z = norm.ppf(1 - alpha / 2)
    tilde_n = n + z**2
    tilde_x = x + (z**2)/2
    tilde_p = tilde_x / tilde_n
    deviation = z * np.sqrt(tilde_p * (1 - tilde_p) / tilde_n)
    lower = tilde_p - deviation
    upper = tilde_p + deviation
    return (lower, upper)


def main():
    n = 20
    alpha = 0.05

    ps = np.arange(0.02, 1, 0.02)
    wald_coverage_probabilities = np.zeros(len(ps))
    wilson_coverage_probabilies = np.zeros(len(ps))
    agresti_coull_coverage_probabilies = np.zeros(len(ps))

    for index, p in enumerate(ps):

        for x in np.arange(n + 1):
            px = binom.pmf(x, n, p)

            wald_ci = wald_confidence_interval(n, x, alpha)
            if wald_ci[0] <= p and p <= wald_ci[1]:
                wald_coverage_probabilities[index] += px

            wilson_ci = wilson_confidence_interval(n, x, alpha)
            if wilson_ci[0] <= p and p <= wilson_ci[1]:
                wilson_coverage_probabilies[index] += px

            agresti_coull_ci = agresti_coull_confidence_interval(n, x, alpha)
            if agresti_coull_ci[0] <= p and p <= agresti_coull_ci[1]:
                agresti_coull_coverage_probabilies[index] += px

    plt.plot(ps, (1 - alpha) * np.ones(len(ps)), linestyle="-", color="black")

    plt.plot(ps, wald_coverage_probabilities, label="wald", linestyle=":", marker="o", color="blue", markerfacecolor='none')

    plt.plot(ps, wilson_coverage_probabilies, label="wilson", linestyle=":", marker="x", color="red")

    plt.plot(ps, agresti_coull_coverage_probabilies, label="agresti-coull", linestyle=":", marker="+", color="green")

    plt.xlabel("$p$")
    plt.ylabel("coverage probability")
    plt.title(f"Coverage Probabilities (n={n} / confidence level={1-alpha})")
    plt.legend(loc="lower center")
    plt.show()


if __name__ == "__main__":
    main()

横軸を真の母比率$p$の値, 縦軸を被覆確率としてプロットするとこのようになる（参考までに, 被覆確率$1-\alpha$のラインも黒の実線で示している）:

グラフから読み取れることとしては以下のとおり:

• Wald interval の被覆確率は, 特に$p=0, 1$の両端に近い部分で急激に低くなる
• Agresti-Coull interval の被覆確率が Wilson intervalの被覆確率以上となる
  • これは先述した「Agresti-Coull interval は Wilson interval を含む」という事実から導かれる

信頼区間の特徴

Wald intervalは$X = 0, n$において信頼区間が1点のみ

• これが成り立つことは, Wald intervalの式に$\hat{p} = 0, 1$を代入してみれば明らか
• 信頼区間が1点のみとなってしまう性質は, 望ましいとは言えない
  • 例えば$n=10, X=0$だった場合, 信頼区間は$0$の1点のみとなる
  • 薬の副作用を調べる際に「10人調べて副作用出た人が0人だったから, 副作用発生率の信頼区間は$0$です」と言われて納得できるだろうか?

Wilson score intervalは0以上1以下の範囲に収まる

• 別途0, 1の境界値の処理を考慮する必要がないのは利点
• これが成り立つことは, 簡単な計算で確かめることができる

（証明はこちら）

• Wilson score intervalの下限が0以上
• Wilson score intervalの上限が1以下

の2点を示せばよい.

まずは前者(Wilson score intervalの下限が0以上)を示す. 示すべき式は

$$\frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} - \frac{z_{\alpha/2}}{\sqrt{n}} \sqrt{\hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n}}}{1 + \frac{z_{\alpha/2}^{2}}{n}} \ge 0 \\ \Leftrightarrow \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \ge \frac{z_{\alpha/2}}{\sqrt{n}} \sqrt{\hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n}}$$

である. そこで, 左辺の2乗 - 右辺の2乗を計算すると

$$\hat{p}^{2} + \frac{z_{\alpha/2}^{2}}{n} \hat{p} + \frac{z_{\alpha/2}^{4}}{4 n^2} - \frac{z_{\alpha/2}^{2}}{n} \hat{p} \left( 1 - \hat{p} \right) - \frac{z_{\alpha/2}^{4}}{4 n^{2}} = \left( 1 + \frac{z_{\alpha/2}^{2}}{n} \right) \hat{p}^{2} \ge 0$$

となり, 左辺$\ge$右辺が示せた. これでWilson score intervalの下限が0以上となることが示せた.

次に後者(Wilson score intervalの上限が1以下)を示す. 示すべき式は

$$\frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} + \frac{z_{\alpha/2}}{\sqrt{n}} \sqrt{\hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n}}}{1 + \frac{z_{\alpha/2}^{2}}{n}} \le 1 \\ \Leftrightarrow \frac{z_{\alpha/2}}{\sqrt{n}} \sqrt{\hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n}} \le 1 - \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n}$$

である. そこで, 右辺の2乗 - 左辺の2乗を計算すると

$$\left( 1 - \hat{p} \right)^{2} + \frac{z_{\alpha/2}^{2}}{n} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{4}}{4 n^2} - \frac{z_{\alpha/2}^{2}}{n} \hat{p} \left( 1 - \hat{p} \right) - \frac{z_{\alpha/2}^{4}}{4 n^{2}} = \left( 1 + \frac{z_{\alpha/2}^{2}}{n} \right) \left( 1 - \hat{p} \right)^{2} \ge 0$$

となり, 左辺$\le$右辺が示せた. これでWilson score intervalの上限が1以下となることが示せた.

（証明終わり）

Agresti–Coull interval は Wilson interval を含む

• ここから「Agresti–Coull interval の被覆確率は Wilson interval の被覆確率以上である」ことが言える
• ざっくりいえば, Agresti–Coull interval の方が Wilson interval より保守的な信頼区間だといえる
• これも, 簡単な計算で確かめることができる

（証明はこちら）

両信頼区間を再掲する:

$$\begin{aligned} CI_{\rm{Wilson}} &= \frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \pm \frac{z_{\alpha/2}}{\sqrt{n}} \sqrt{\hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n}}}{1 + \frac{z_{\alpha/2}^{2}}{n}}, \\ CI_{\rm{Agresti–Coull}} &= \frac{\hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \pm \frac{z_{\alpha/2}}{\sqrt{n + z_{\alpha/2}^{2}}} \sqrt{\left( \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \right) \left( 1 - \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \right)}}{1 + \frac{z_{\alpha/2}^{2}}{n}}. \end{aligned}$$

両区間とも中央の値は同じため, 幅に相当する部分の大小関係によって信頼区間の包含関係が決まる. そこで

$$\begin{aligned} d_{\rm{Wilson}} &= \frac{z_{\alpha/2}}{\sqrt{n}} \sqrt{\hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n}}, \\ d_{\rm{Agresti–Coull}} &= \frac{z_{\alpha/2}}{\sqrt{n + z_{\alpha/2}^{2}}} \sqrt{\left( \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \right) \left( 1 - \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \right)} \end{aligned}$$

と置き, $d_{\rm{Agresti–Coull}} \ge d_{\rm{Wilson}}$を示すことで$CI_{\rm{Agresti–Coull}} \supseteq CI_{\rm{Wilson}}$を示す.

2乗の差をとると

$$\begin{aligned} & d_{\rm{Agresti–Coull}}^{2} - d_{\rm{Wilson}}^{2} \\ = & \frac{z_{\alpha/2}^{2}}{n + z_{\alpha/2}^{2}} \left( \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \right) \left( 1 - \hat{p} + \frac{z_{\alpha/2}^{2}}{2 n} \right) - \frac{z_{\alpha/2}^{2}}{n} \left( \hat{p} \left( 1 - \hat{p} \right) + \frac{z_{\alpha/2}^{2}}{4 n} \right) \\ = & \frac{z_{\alpha/2}^{4}}{n \left( n + z_{\alpha/2}^{2} \right)} \left( \frac{1}{4} - \hat{p} \left( 1 - \hat{p} \right) \right) \\ \ge & 0 \end{aligned}$$

となる($\hat{p} = \frac{1}{2}$で等号成立). 従って, $d_{\rm{Agresti–Coull}} \ge d_{\rm{Wilson}}$ すなわち $CI_{\rm{Agresti–Coull}} \supseteq CI_{\rm{Wilson}}$ が示された.

（証明終わり）

Wald interval の被覆確率は, 特に$p=0, 1$の両端に近い部分で急激に低くなる

Brown et al.(2002)のsection 2で議論があるため, 詳細はそちらを参照. 論文中では, 本記事で紹介したもの以外のものも含めた複数の信頼区間の被覆確率をEdgeworth展開し, その結果に基づいた比較が行われている.

個人的な所感

• あえて Wald interval を推す理由はなさそう
  • 理論的な正当化はできるが, 統計初学者にはハードルが高い
  • 被覆確率の挙動も望ましくない（特に$p=0, 1$周辺）
  • とはいえ計算式は非常にシンプルなので, 統計学の試験など, 計算資源が乏しい中短時間で計算をしないといけないような限定的な場面であれば, 使う意義はありそう
• 実際に信頼区間を計算する際は, 以下のように使い分けをするとよさそう
  • 基本的には Wilson interval を利用
    • 3つの中で式は一番複雑に見えるが, 導出の過程は自然に感じられた
    • 被覆確率を見てもさほど大きな問題はなく, 信頼区間が0以上1以下の範囲に収まる点も良い
  • 簡便さを重視する場合は Agresti-Coull interval を利用
    • Wald interval のような式の簡便さと被覆確率の安定度合いを両方兼ね備えている
    • ただ, 信頼区間が0以上1以下の範囲に収まるとは限らず, 多少のケアは必要
• 母比率の信頼区間の使い分けについては Brown et al.(2001) の議論が参考になる

参考文献

• Brown, L. D., Cai, T. T., and DasGupta, A. (2001). Interval Estimation for a Binomial Proportion. Statistical Science, 16(2), 101-133.
• Brown, L. D., Cai, T. T., and DasGupta, A. (2002). Confidence Intervals for a binomial proportion and asymptotic expansions. The Annals of Statistics, 30(1), 160–201