バートレット検定

前回の記事の中で、SciPyにはF検定の関数が実装されていないという話を書きました。
参考: PythonでF検定を実装する

しかし、その代わりにSciPyには分散が等しいことを検定する別の方法が実装されています。
その一つが、バートレット検定です。
関数はこれです。
scipy.stats.bartlett

複数(２個だけでなく3個以上も可能)のサンプルを渡してあげると、統計量とp値を計算してくれます。

今回の記事では、使い方の前にバートレット検定の統計量を紹介します。
参考にしたのは、東京大学出版会の自然科学の統計学 (青い本) の 94ページです。

サンプルの数を$a$とし、それぞれのサンプルのサイズを$n_i \ \ (j=1,\dots, a)$とします。
また、サンプルサイズの合計を$n=\sum_{j}n_j$と置いておきます。
各サンプルの不偏分散を
$$
V_i = \sum_{j} \frac{(y_{ij}-\bar{y_i})^2}{n_i-1}
$$
として、さらにこれらを併合したものを
$$
V_e = \sum_i \frac{(n_i-1)V_i}{n-a} = \sum_i \sum_j \frac{(y_{ij}-\bar{y_i})^2}{n-a}
$$
とします。
さらに、
$$
B = (n-a)\log{V_e} -\sum_i (n_i-1)\log{V_i}
$$
と置いて、
$$
B’ = \frac{B}{
1+\frac{1}{3(a-1)}\{\sum_i \frac{1}{n_i-1} – \frac{1}{n-a}\}
}
$$
と補正すると、この$B’$は自由度$a-1$の$\chi^2$分布に近似的にしたがいます。
それを利用して、検定を行うのがバートレット検定です。

数式に沿ってNumPyで実装してみたのが次の関数です。ランダム生成した３サンプルで試してみました。

from scipy.stats import norm
from scipy.stats import chi2
import numpy as np

# スクラッチで実装したバーレット検定量の算出関数
def bartlett_value(*args):
    # サンプル数
    a = len(args)
    # 各サンプルの　サンプルサイズと不偏分散
    ni = np.zeros(a)
    Vi = np.zeros(a)
    for j in range(a):
        ni[j] = len(args[j])
        Vi[j] = np.var(args[j], ddof=1)

    # サンプルサイズの合計
    N = np.sum(ni)
    # Veの計算
    Ve = np.sum((ni-1)*Vi) / (N-a)

    # 統計量の計算
    B = (N-a) * np.log(Ve) - np.sum((ni-1)*np.log(Vi))
    C = 1+1/(3*(a-1))*(np.sum(1/(ni-1))-1/(N-a))

    return B/C


# サンプルを３種類生成する
data_a = norm(loc=3, scale=5).rvs(20)
data_b = norm(loc=7, scale=8).rvs(30)
data_c = norm(loc=-5, scale=5).rvs(15)

# 統計量を算出
b_value = bartlett_value(data_a, data_b, data_c)

# p値を算出
p_value = chi2(3-1).sf(b_value)

print("B値:", b_value)
# B値: 8.26319721031961
print("p値:", p_value)
# p値: 0.016057189188779606

$\alpha=0.05$とすると帰無仮説が棄却され、分散が等しくないことが検出されていますね。

もちろん、SciPyに関数は用意されているので、普段はこの様にスクラッチで実装する必要はなく、
呼び出すだけで良いです。結果が一致することを見ておきます。

from scipy.stats import bartlett


b_value, p_value = bartlett(data_a, data_b, data_c)
print("B値:", b_value)
# B値: 8.26319721031961
print("p値:", p_value)
# p値: 0.016057189188779606

バッチリ一致しました。