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

これは何か

• 母比率の信頼区間について調べたことをまとめたもの
• 信頼区間構築のやり方はいくつもあるが, 本記事では正規分布への近似を利用しないものについてまとめた
  • 具体的には, 以下の2つの信頼区間について整理した
    • Clopper-Pearson interval
    • Jeffreys interval^[1]
• 正規分布への近似を利用したやり方についてはこちらの記事でまとめた

記号の定義

• Binom$(n, p)$: 2項分布
• Beta$(a, b)$: ベータ分布
  • $\Beta \left(a, b \right)$: ベータ関数
  • $\Beta \left(\alpha; a, b \right)$: ベータ分布B$(a, b)$の下側$100\alpha$%点

各信頼区間の概要

以後, $X \sim$ Binom$(n, p)$とする.

Clopper-Pearson interval

• 2項分布に基づく片側検定を2つ考えることにより導かれる信頼区間

$$CI_{{\rm Clopper-Pearson}} = \left[ \Beta \left(\frac{\alpha}{2}; X, n − X + 1 \right), \Beta \left(1-\frac{\alpha}{2}; X + 1, n − X \right) \right].$$

ただし, $X=0$の時の下限は0, $X=n$の時の上限は1とする.

Clopper-Pearson intervalの導出

信頼区間の上限・下限をそれぞれ導出していく.

信頼区間の下限を導くにあたり, 有意水準$100 \times \frac{\alpha}{2}$%の片側検定$H_{0}: p = p_{0}$ v.s. $H_{1}: p > p_{0}$を考える. 有意水準を2で割っているのは, 後ほどもう1つの片側検定を考えるためである. 帰無仮説が棄却されない$p_{0}$の範囲は

$$\sum_{k=X}^{n} \binom{n}{k} p_{0}^{k} \left( 1 - p_{0} \right)^{n - k} \gt \frac{\alpha}{2}$$

となる.

$X=0$ の場合, 左辺は1となり, $p_{0}$の値によらず不等式が成立してしまう. この場合は$p_{0}$の下限値を0と定義する. 下限の導出においては, 以後$X \gt 0$を仮定する. 左辺はベータ分布Beta$(X, n - X + 1)$の下側確率で表現されるため, この事実を使うと解を比較的シンプルな形で書き下すことができる. 左辺がベータ分布の下側確率で表現されることは, 部分積分を繰り返し適用することで容易に示せる:

$$\begin{aligned} & \frac{1}{\Beta \left( X, n - X + 1 \right)} \int_{0}^{p_{0}} p^{X - 1} \left( 1 - p \right)^{n - X} {\rm d}p \\ =& \int_{0}^{p_{0}} \frac{\Gamma \left( n + 1 \right)}{\Gamma \left( X \right) \Gamma \left( n - X + 1 \right)} p^{X - 1} \left( 1 - p \right)^{n - X} {\rm d}p \\ =& \int_{0}^{p_{0}} \frac{n!}{\left( X - 1 \right)! \left( n - X \right)!} p^{X - 1} \left( 1 - p \right)^{n - X} {\rm d}p \\ =& \frac{n!}{X! \left( n - X \right)!} p_{0}^{X} \left( 1 - p_{0} \right)^{n - X} + \int_{0}^{p_{0}} \frac{n!}{X! \left( n - X - 1 \right)!} p^{X} \left( 1 - p \right)^{n - X - 1} {\rm d}p \\ =& \cdots \\ =& \sum_{k=X}^{n-1} \frac{n!}{k! \left( n - k \right)!} p_{0}^{k} \left( 1 - p_{0} \right)^{n - k} + \int_{0}^{p_{0}} \frac{n!}{\left( n - 1 \right)!} p^{n-1} {\rm d}p \\ =& \sum_{k=X}^{n} \frac{n!}{k! \left( n - k \right)!} p_{0}^{k} \left( 1 - p_{0} \right)^{n - k}. \end{aligned}$$

つまり元の不等式を言い換えると, 「ベータ分布の下側確率が$\frac{\alpha}{2}$より大きい」ということになる. ゆえに

$$p_{0} \gt \Beta \left(\frac{\alpha}{2}; X, n − X + 1 \right)$$

が不等式の解となる. これで求める信頼区間の下限の値が求められた.

上限の値についても同様に考える. 有意水準$100 \times \frac{\alpha}{2}$%の片側検定$H_{0}: p = p_{0}$ v.s. $H_{1}: p < p_{0}$を考える（対立仮説の不等号の向きが先ほどと逆であることに注意）. すると, 帰無仮説が棄却されない$p_{0}$の範囲は

$$\sum_{k=0}^{X} \binom{n}{k} p_{0}^{k} \left( 1 - p_{0} \right)^{n - k} \gt \frac{\alpha}{2}$$

となる.

$X=n$ の場合, 左辺は1となり, $p_{0}$の値によらず不等式が成立してしまう. この場合は$p_{0}$の上限値を1と定義する. 上限の導出においては, 以後$X \lt n$を仮定する. この不等式の左辺は, ベータ分布Beta$(X + 1, n - X)$の上側確率で表現される（先ほどの式を利用）:

$$\begin{aligned} & \sum_{k=0}^{X} \frac{n!}{k! \left( n - k \right)!} p_{0}^{k} \left( 1 - p_{0} \right)^{n - k} \\ &= 1 - \sum_{k=X + 1}^{n} \frac{n!}{k! \left( n - k \right)!} p_{0}^{k} \left( 1 - p_{0} \right)^{n - k} \\ &= 1 - \frac{1}{\Beta \left( X + 1, n - X \right)} \int_{0}^{p_{0}} p^{X} \left( 1 - p \right)^{n - X - 1} {\rm d}p \\ &= \frac{1}{\Beta \left( X + 1, n - X \right)} \int_{p_{0}}^{1} p^{X} \left( 1 - p \right)^{n - X - 1} {\rm d}p. \end{aligned}$$

つまり元の不等式を言い換えると, 「ベータ分布の上側確率が$\frac{\alpha}{2}$より大きい」, 「ベータ分布の下側確率が$1 - \frac{\alpha}{2}$以下」ということになる. ゆえに

$$p_{0} \le \Beta \left( 1 - \frac{\alpha}{2}; X + 1, n − X \right)$$

が不等式の解となる. これで求める信頼区間の上限の値が求められた.

以上をまとめて, 信頼水準$100(1-\alpha)$%のClopper-Pearson信頼区間は

$$\left[ \Beta \left(\frac{\alpha}{2}; X, n − X + 1 \right), \Beta \left(1-\frac{\alpha}{2}; X + 1, n − X \right) \right]$$

となる. ただし, $X=0$の時の下限は0, $X=n$の時の上限は1とする.

Jeffreys interval

• ベイズ統計の枠組みで導出される信用区間
• 2項分布モデルにJeffreys priorを適用して得られる事後分布から計算される

$$CI_{{\rm Jeffreys}} = \left[ \Beta \left(\frac{\alpha}{2}; X + \frac{1}{2}, n − X + \frac{1}{2} \right), \Beta \left(1-\frac{\alpha}{2}; X + \frac{1}{2}, n − X + \frac{1}{2} \right) \right].$$

Jeffreys intervalの導出

2項分布モデルBinom$(n, p)$に対し, 事前分布としてJeffreys priorを採用した際の事後分布から信用区間を算出する.

Binom$(n, p)$のJeffreys priorを計算すると,

$$p_{{\rm Jeffreys}} \left( p \right) \propto p^{-\frac{1}{2}} \left( 1 - p \right)^{-\frac{1}{2}},$$

つまりベータ分布Beta$\left( \frac{1}{2}, \frac{1}{2} \right)$となる.

Jeffreys priorの導出

Jeffreys priorは, Fisher情報量を$I\left(p\right)$として

$$p_{{\rm Jeffreys}} \left( p \right) \propto \left| I\left(p\right) \right|^{-\frac{1}{2}}$$

で定義される. そのため, 2項分布のFisher情報量を計算する. 2項分布のスコア関数（対数尤度のパラメータに関する1階偏微分）は

$$\begin{aligned} \frac{\partial}{\partial p} \log p\left(X | p\right) &= \frac{X}{p} - \frac{n - X}{1 - p} \\ &= \frac{X - np}{p (1 - p)} \end{aligned}$$

なので, Fisher情報量は

$$\begin{aligned} I\left(p\right) &= E\left[ \left( \frac{\partial}{\partial p} \log p\left(X | p\right) \right)^{2} \right] \\ &= E\left[ \left( \frac{X - np}{p (1 - p)} \right)^{2} \right] \\ &= \frac{V \left[ X \right]}{p^2 (1-p)^2} \\ &= \frac{n}{p (1 - p)} \end{aligned}$$

となる.

ゆえにJeffreys priorは

$$p_{{\rm Jeffreys}} \left( p \right) \propto \left| \frac{n}{p \left(1 - p \right)} \right|^{-\frac{1}{2}} \propto p^{-\frac{1}{2}} \left(1 - p \right)^{-\frac{1}{2}}$$

となる. （Jeffreys priorの導出終わり）

この事前分布を用いると, $p$の事後分布は

$$\begin{aligned} p \left( p | X \right) & \propto p \left( X | p \right) \times p_{{\rm Jeffreys}} \left( p \right) \\ & \propto p^{X} \left( 1 - p \right)^{n-X} p^{-\frac{1}{2}} \left( 1 - p \right)^{-\frac{1}{2}} \\ & = \propto p^{X-\frac{1}{2}} \left( 1 - p \right)^{n-X-\frac{1}{2}} \end{aligned}$$

つまりベータ分布Beta$\left( X + \frac{1}{2}, n - X + \frac{1}{2} \right)$となる.

あとは, 両端から$100 \times \frac{\alpha}{2}$%の距離にある2点を上限・下限とすれば, Jeffreys interval

$$\left[ \Beta \left(\frac{\alpha}{2}; X + \frac{1}{2}, n − X + \frac{1}{2} \right), \Beta \left(1-\frac{\alpha}{2}; X + \frac{1}{2}, n − X + \frac{1}{2} \right) \right]$$

が求まる.

信頼区間の観察

参考までに, こちらの記事で紹介したAgresti–Coull intervalと合わせて可視化を行う.

信頼区間の数値例

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

（可視化コード）

confidence_intervals.py

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


def clopper_pearson_confidence_interval(n: int, x: int, alpha: float = 0.05) -> tuple[float, float]:
    lower, upper = 0, 1
    if x > 0:
        lower = beta.ppf(alpha / 2, x, n - x + 1)
    if x < n:
        upper = beta.ppf(1 - alpha / 2, x + 1, n - x)
    return (lower, upper) 


def jeffreys_confidence_interval(n: int, x: int, alpha: float = 0.05) -> tuple[float, float]:
    lower = beta.ppf(alpha / 2, x + 1/2, n - x + 1/2)
    upper = beta.ppf(1 - alpha / 2, x + 1/2, n - x + 1/2)
    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)
    clopper_pearson_lower, clopper_pearson_upper = np.zeros(n + 1), np.zeros(n + 1)
    jeffreys_lower, jeffreys_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:
        clopper_pearson_ci = clopper_pearson_confidence_interval(n, x, alpha)
        clopper_pearson_lower[x], clopper_pearson_upper[x] = clopper_pearson_ci[0], clopper_pearson_ci[1]

        jeffreys_ci = jeffreys_confidence_interval(n, x, alpha)
        jeffreys_lower[x], jeffreys_upper[x] = jeffreys_ci[0], jeffreys_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, clopper_pearson_lower, label="clopper-pearson", linestyle=":", marker="o", color="blueviolet", markerfacecolor='none')
    plt.plot(xs, clopper_pearson_upper, linestyle=":", marker="o", color="blueviolet", markerfacecolor='none')

    plt.plot(xs, jeffreys_lower, label="jeffreys", linestyle=":", marker="x", color="orange")
    plt.plot(xs, jeffreys_upper, linestyle=":", marker="x", color="orange")

    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$の値, 縦軸を母比率の信頼区間の上限・下限の値としてプロットしている:

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

• Clopper-Pearson intervalは他のintervalsと比べて区間が広い
  • この特徴は後述する被覆確率にも影響してくる
• Clopper-Pearson interval / Jeffreys intervalは0以上1以下の範囲に収まる
  • こちらの記事でも触れたが, Agresti-Coull intervalなどでは0以上1以下の範囲からはみ出ることもある

信頼区間の被覆確率

被覆確率を「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}.$$

他の母比率の信頼区間の解説記事を見ると, 被覆確率をモンテカルロシミュレーションで計算して比較しているケースが多い印象があるが, 被覆確率については上記のようにシミュレーションなしで計算が可能である.

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

（可視化コード）

coverage_probabilities.py

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


def clopper_pearson_confidence_interval(n: int, x: int, alpha: float = 0.05) -> tuple[float, float]:
    lower, upper = 0, 1
    if x > 0:
        lower = beta.ppf(alpha / 2, x, n - x + 1)
    if x < n:
        upper = beta.ppf(1 - alpha / 2, x + 1, n - x)
    return (lower, upper) 


def jeffreys_confidence_interval(n: int, x: int, alpha: float = 0.05) -> tuple[float, float]:
    lower = beta.ppf(alpha / 2, x + 1/2, n - x + 1/2)
    upper = beta.ppf(1 - alpha / 2, x + 1/2, n - x + 1/2)
    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)
    clopper_pearson_coverage_probabilities = np.zeros(len(ps))
    jeffreys_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)

            clopper_pearson_ci = clopper_pearson_confidence_interval(n, x, alpha)
            if clopper_pearson_ci[0] <= p and p <= clopper_pearson_ci[1]:
                clopper_pearson_coverage_probabilities[index] += px

            jeffreys_ci = jeffreys_confidence_interval(n, x, alpha)
            if jeffreys_ci[0] <= p and p <= jeffreys_ci[1]:
                jeffreys_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, clopper_pearson_coverage_probabilities, label="clopper-pearson", linestyle=":", marker="o", color="blueviolet", markerfacecolor='none')

    plt.plot(ps, jeffreys_coverage_probabilies, label="jeffreys", linestyle=":", marker="x", color="orange")

    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$のラインも黒の実線で示している）:

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

• Jeffreys intervalは, 被覆確率が想定している被覆確率のラインの周辺に分布している
• Clopper-Pearson intervalは, 想定している被覆確率のラインと比べてそれなりに高い被覆確率を保っている

個人的な所感

• Clopper-Pearson intervalの被覆確率の高さは活用シーンを選びそう
  • Brown et al.(2001)のSection 4.2では wastefully conservativeと評されている
  • とはいえ保守的に寄った意思決定をしたいシーンでは活用できるはず

参考文献

• 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

脚注

1. これは信用区間に分類されるものではあるが, 本記事では特に区別はしない