2項検定

Rを使った検定と区間推定もご参照ください。

2項分布

歪んだ硬貨があって、投げると確率 $\theta$ で表（おもて）が出て、確率 $1 - \theta$ で裏が出るとします。この硬貨を $n$ 回投げて、表が $r$ 回出る確率を求めてください。

これが2項分布（binomial distribution）の問題です。答えは

\[ _nC_r \theta^r (1-\theta)^{n-r} \]

です。$_nC_r$ は $n$ 個から $r$ 個を選ぶ組合せ（combination）の数です。階乗（!）を使えば

\[ _nC_r = \frac{n!}{r!(n-r)!} \]

と表せます。

Python では、階乗は math.factorial() という関数ですので、例えば $_{10}C_3$ は次のようにして求められます。

import math

math.factorial(10) // (math.factorial(3) * math.factorial(7))

120

ここで // は整数の除算を意味します。

上の計算は math.factorial(10) / (math.factorial(3) * math.factorial(7)) でも同じですが、// を使えば結果は 120 なのに対し、/ を使えば結果は 120.0 になります。つまり / は分子・分母が整数でも浮動小数点数に直してから除算します。Python はどんな大きな整数でも正確に計算できるのに対し、浮動小数点数は CPU で決まった精度でしか計算できません。例えば math.factorial(200) // 10 は正確に計算できますが math.factorial(200) / 10 は分子が大きすぎて浮動小数点数に直せないのでエラーになります。

Python 3.8 以降では math.comb(10, 3) でも求められます。こちらのほうが高速です。ただ、Coogle Colaboratory がまだ Python 3.7 のようです。3.7 以前では次のようにして定義すると便利です。

def comb(n, r):
    return math.factorial(n) // (math.factorial(r) * math.factorial(n - r))

投げると確率 0.4 で表が出る硬貨を10回投げて表が3回出る確率 $_{10}C_3 \cdot 0.4^3 \cdot 0.6^7$ は

math.factorial(10) // (math.factorial(3) * math.factorial(7)) * 0.4**3 * 0.6**7

0.21499084799999998

です（正確な値は 0.214990848 です）。

同じ値は SciPy ライブラリを使っても求められます。SciPy ライブラリは非常に大きいので、統計の部分 stats だけインポートします。

from scipy import stats

stats.binom.pmf(3, 10, 0.4)

0.21499084799999982

stats.binom.pmf() は binomial probability mass function の意味です。

投げると確率 0.5 で表が出る硬貨を10回投げて表が0〜10回出る確率を全部出力するには次のようにします。

stats.binom.pmf([0,1,2,3,4,5,6,7,8,9,10], 10, 0.5)

array([0.00097656, 0.00976563, 0.04394531, 0.1171875 , 0.20507812,
       0.24609375, 0.20507812, 0.1171875 , 0.04394531, 0.00976563,
       0.00097656])

0 以上 11 未満の整数 [0,1,2,3,4,5,6,7,8,9,10] は range(11) とも書けます。

これをグラフにするには次のように打ち込みます。

import matplotlib.pyplot as plt

plt.bar(range(11), stats.binom.pmf(range(11), 10, 0.5))
グラフ

確率なので、全部加えると 1 になります。

sum(stats.binom.pmf(range(11), 10, 0.5))

1.0000000000000004

ちなみに、stats.binom.pmf() の累積を求める stats.binom.cdf() という関数もあります。cdf は cumulative distribution function の意味です。

stats.binom.cdf(range(11), 10, 0.5)

array([9.76562500e-04, 1.07421875e-02, 5.46875000e-02, 1.71875000e-01,
       3.76953125e-01, 6.23046875e-01, 8.28125000e-01, 9.45312500e-01,
       9.89257812e-01, 9.99023438e-01, 1.00000000e+00])

統計的検定の考え方

授業ではJavaScriptのプログラムも使って説明します。

硬貨を10回投げて表が2回しか出ませんでした。この硬貨は歪んでいるでしょうか。

ある母集団（ここでは例えばある高校の全生徒のように人数の多い集団を考えます）から10人をランダムに選んで聞いたところ、賛成2人、反対8人でした。母集団全体でも反対のほうが多いと言えるでしょうか。

これらの問いについて考えるために、仮に硬貨は歪んでいない（あるいは母集団全体では賛否が等しい）というモデル（帰無仮説、null hypothesis）を設定します。そして、この帰無仮説が正しかった場合に、実際に観測された事象、およびそれより珍しい事象の確率の合計を求めます。この場合は、表が0〜2回・8〜10回出る確率の合計です。この確率の合計を $p$ 値（ピーち、$p$-value）といいます。

$p$ 値が非常に小さければ、実際に起きた事象はこのモデルでは説明しにくいので、たぶん硬貨は歪んでいる（あるいは賛否は等しくない）と推測します。$p$ 値が大きければ、これだけのデータでは何も言えないので、判断を保留します（硬貨が歪んでない・賛否が等しいと推測するわけではありません）。

$p$ 値が大きいか小さいかの境界（有意水準）を仮に 0.05 として、$p \leqq 0.05$（または $p < 0.05$）であれば帰無仮説からの外れが「統計的に有意」（statistically significant）である、あるいは「帰無仮説は棄却（reject）される」ということがあります。0.05 という値に特に意味はありませんが、伝統的によく使われています（物理学では通常もっともっと厳しい条件を課します）。

さて、硬貨投げで表も裏も 1/2 の確率で出るという帰無仮説では、表が $r$ 回出る確率は $_{10}C_r (1/2)^{10}$ です。表が0、1、2、8、9、10回出る確率はそれぞれ

stats.binom.pmf([0,1,2,8,9,10], 10, 0.5)

array([0.00097656, 0.00976563, 0.04394531, 0.04394531, 0.00976563,
       0.00097656])

で、この合計、すなわち $p$ 値は

sum(stats.binom.pmf([0,1,2,8,9,10], 10, 0.5))

0.1093750000000001

になります。同じことが

stats.binom.cdf(2, 10, 0.5) * 2

0.109375

でも求められます。また、次のようにすれば、このような計算（2項検定）が簡単にできます。

stats.binomtest(2, 10, 0.5).pvalue

0.109375

Python では SciPy に，$p$ 値を返すだけの古い scipy.stats.binom_test と，2021-06-20 の SciPy 1.7.0 で新しく追加された scipy.stats.binomtest とがあります。計算量は，前者は $O(n)$，後者は2分探索を使うので $O(\log n)$ だそうです。古い方は次のように使います。

  
stats.binom_test(2, 10, 0.5)

binom_test() はSciPy 1.10.0の段階で「deprecated」で、1.12.0で削除される予定です。

つまり、表が2回出るという事象の $p$ 値は 0.1094 ほどです。したがって、有意水準を 0.05 とすれば、表裏の差は統計的に有意ではありませんし、アンケートであればこんなに少人数の結果から「賛成が少ない」という結論を導いてはいけないということになります。

表が1回（賛成が1人）なら、$p$ 値は 0.02 ほどになり、水準 0.05 で有意になります。

これが、フィッシャー（R. A. Fisher、1890〜1962年）が「有意性の検定」（tests of significance, significance tests）と呼んだ方法の考え方です。

ここでは、実際に生じた事象またはそれより珍しい事象の生じる確率の合計として $p$ 値を定義しました。これは2項分布のような山の形の分布では両側の部分の確率になりますので、両側検定といいます。一方、片側だけを考える片側検定もあります。例えば表が2回出た場合の（左側だけの）片側 $p$ 値は、表が0回・1回・2回出る確率の合計です。

ICH（医薬品規制調和国際会議）の定めるICH E9「臨床試験のための統計的原則」というガイドラインでは、「原則として片側仮説を検証する場合は2.5%、両側仮説の場合は5%とすること」とされています。

区間推定

授業ではJavaScriptのプログラムも使って説明します。

なお、高校「数学B」では、2項分布 $B(n, \theta)$ の平均が $n\theta$、分散が $n\theta(1-\theta)$ であることを使って、同じ平均と分散を持つ正規分布で近似して扱うことが多いと思います（中心極限定理から $n \to \infty$ で正規分布になります）。

10回投げて2回表が出たとき、表が出る確率の95%信頼区間は次のようにして求めます。これは正規分布近似ではなく Clopper-Pearson の正確な方法を使っています。

stats.binomtest(2, 10).proportion_ci()

ConfidenceInterval(low=0.025210726326833497, high=0.5560954623076415)

[0.025, 0.556] であることがわかります。

下限・上限を別々に求めるには次のようにします。

ci = stats.binomtest(2, 10).proportion_ci()
print(ci.low, ci.high)

0.025210726326833497 0.5560954623076415

あるいは次のようにもできます。

lo, hi = stats.binomtest(2, 10).proportion_ci()
print(lo, hi)

信頼区間は statsmodels.stats.proportion.proportion_confint でも求められます：

from statsmodels.stats.proportion import proportion_confint

proportion_confint(2, 10, method="beta")

(0.02521072632683336, 0.5560954623076415)

Rの binom.test() をPythonで実現

Rには binom.test() という関数があります：

binom.test(x, n, p = 0.5,
           alternative = c("two.sided", "less", "greater"),
           conf.level = 0.95)

実行例：

> binom.test(2, 10)

	Exact binomial test

data:  2 and 10
number of successes = 2, number of trials = 10, p-value = 0.1094
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.02521073 0.55609546
sample estimates:
probability of success 
                   0.2

これとほぼ同じことをするPythonの関数を作ってみました。

from scipy import stats

def binom_test(x, n, p=0.5, alternative="two-sided", conf_level=0.95):
    r = stats.binomtest(x, n, p, alternative)
    print()
    print("	Exact binomial test")
    print()
    print(f"data:  {x} and {n}")
    print(f"number of successes = {x}, number of trials = {n}, p-value = {r.pvalue:.4g}")
    if alternative == "two-sided":
        s = "not equal to"
    else:
        s = alternative + " than"
    print(f"alternative hypothesis: true probability of success is {s} {p}")
    print(f"{conf_level * 100:g} percent confidence interval:")
    ci = r.proportion_ci(confidence_level=conf_level)
    print(f" {ci.low:.8f} {ci.high:.8f}")
    print("sample estimates:")
    print("probability of success ")
    print(f"                   {r.proportion_estimate:.4g} ")

実行例：

binom_test(2, 10)


	Exact binomial test

data:  2 and 10
number of successes = 2, number of trials = 10, p-value = 0.1094
alternative hypothesis: true probability of success is not equal to 0.5
95.0 percent confidence interval:
 0.02521073 0.55609546
sample estimates:
probability of success 
                   0.2