はじめに
Pythonでt検定したい場合、 scipy.stats
のttest_ind
を使うことが多いと思います。
基本的にはttest_ind
を使用すれば良いと思うのですが、この関数は全サンプルが必要になります。小規模のデータサイズなら何も考えず突っ込めば良いのですが、これが10~100万となったらしんどいですよね。というか、t検定では要約統計量(平均値や不偏分散など)のみが必要なので、わざわざ全サンプル用意するのは無駄な気がしてしまいます。
というわけで、準備するのが要約統計量だけで済むよう、t検定をPythonでスクラッチ実装してみます。scipyのttest_ind
と比較することで、妥当性を検証しようと思います。
2022.01.20追記
scipyで ttest_ind_from_stats という関数が実装されており、要約統計量を引数に渡すことでt検定の結果を得られます。
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind_from_stats.html
そもそもt検定とは?
母平均に対する検定。「グループA,Bの平均値の差」を検定する。
例えば、「A/Bテストを行った結果、Aグループの購買金額の平均は1000円、Bグループの購買金額の平均は1050円だった。この50円の差は偶然の産物なのか?どうなんだ???」ってときに使うイメージです。(違ったらツッコミください)
今回は対応のない2群のt検定を実装します。対応のある2群のt検定、1群のt検定は扱いません。
notebook
githubにコードをアップしています。
https://github.com/hatena-hanata/ml_stats_study/blob/master/stats/01_t_test.ipynb
実装
1. データセット
グループA:平均50、標準偏差10の正規分布から1000個ランダムサンプリングした標本
グループB:平均48、標準偏差12の正規分布から1200個ランダムサンプリングした標本 とします。
今回は正規分布からのランダムサンプリングとしましたが、標本数が大きければ、中心極限定理により標本平均の分布は正規分布に近付くので、標本の分布は正規分布でなくても大丈夫です(間違ってたらご指摘ください)。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
m = 1000
n = 1200
a = np.random.normal(loc=50, scale=10, size=m)
b = np.random.normal(loc=48, scale=12, size=n)
plt.hist(a, bins=100, alpha=0.3, histtype='stepfilled', color='r', label='A')
plt.hist(b, bins=100, alpha=0.3, histtype='stepfilled', color='b', label='B')
plt.legend()
Aの平均値は50.19、Bの平均値は48.58となりました。このAとBの平均値の差が、統計的に有意か調べるのがt検定です。
2. Studentのt検定(等分散を仮定)
まずはじめに、等分散を仮定したStudentのt検定を実施します。
2-1. 各種統計量を求める
# a 平均と不偏分散
a_mean = a.mean()
a_var = a.var(ddof=1)
# b 平均と不偏分散
b_mean = b.mean()
b_var = b.var(ddof=1)
# poolした分散
s_2 = ((m - 1) * a_var + (n - 1) * b_var) / (m + n - 2)
# 自由度
dof = m + n - 2
2-2. t値、p値を求める
from scipy.stats import t
t_value = (a_mean - b_mean) / (np.sqrt(s_2 * (1/m + (1/n))))
dof = m + n - 2
p_value = min(t.cdf(t_value, dof), 1-t.cdf(t_value, dof)) * 2
print(f't値:{t_value} p値:{p_value}')
### output
### t値:3.4287178596465306 p値:0.0006176178409258437
2-3. Scipyのt検定
scipy.statsのttest_ind
を使用します。equal_var=Trueとすることで、Studentのt検定になります。
stats.ttest_ind(a, b)
### output
### Ttest_indResult(statistic=3.4287178596465306, pvalue=0.0006176178409258371)
自分で実装したものと、おおよそ値が一致しています!
3. Weltchのt検定(等分散を仮定しない)
今回のテーマとは関係ないですが、こんな記事を見つけました。
https://biolab.sakura.ne.jp/welch-test.html
等分散かどうかに関わらずWeltchのt検定を使用するのが良さそうです。実装しましょう。
3-1. 各種統計量を求める
Studentのt検定とWeltchのt検定の差異は、分散と自由度が異なることに由来します。
# a 平均と不偏分散
a_mean = a.mean()
a_var = a.var(ddof=1)
# b 平均と不偏分散
b_mean = b.mean()
b_var = b.var(ddof=1)
# 使用する分散
s_2 = a_var / m + b_var / n
# 自由度
dof = int((a_var / m + b_var / n)**2 / (((a_var/m)**2 / (m-1)) + ((b_var/n)**2 / (n-1))))
自由度は少し複雑ですね。
3-2. t値、p値を求める
from scipy.stats import t
t_value = (a_mean - b_mean) / np.sqrt(s_2)
p_value = min(t.cdf(t_value, dof), 1-t.cdf(t_value, dof)) * 2
print(f't値:{t_value} p値:{p_value}')
### ouput
### t値:3.490154519416705 p値:0.000492256894277121
3-3. Scipyのt検定
ttest_ind
にequal_var=False
とすることで、Weltchのt検定になります。
stats.ttest_ind(a, b, equal_var=False)
### output
### Ttest_indResult(statistic=3.490154519416705, pvalue=0.0004922547773957867)
値がおおよそ一致していることが分かります。
関数化してみた
以上の処理を1つの関数にまとめてしまいます。
from scipy.stats import t
import numpy as np
def t_test(a_mean, a_var, a_n, b_mean, b_var, b_n, equal_var=True):
# Studentのt検定
if equal_var:
s_2 = ((a_n - 1) * a_var + (b_n - 1) * b_var) / (a_n + b_n - 2)
dof = a_n + b_n - 2
# Weltchのt検定
else:
s_2 = a_var / a_n + b_var / b_n
dof = int((a_var / a_n + b_var / b_n)**2 / (((a_var/ a_n)**2 / (a_n-1)) + ((b_var/b_n)**2 / (b_n-1))))
t_value = (a_mean - b_mean) / np.sqrt(s_2)
p_value = min(t.cdf(t_value, dof), 1-t.cdf(t_value, dof)) * 2
print(f'TtestResult(statistic={t_value}, pvalue={p_value})')
引数が多いので、使い勝手は正直微妙ですね^^;
まとめ
上の実装の過程を見てもらえれば分かる通り、平均値、不偏分散、サンプルサイズが分かればt検定を実施することができます。
DBの大規模データに対してSQLで要約統計量を計算→Pythonでt検定実施
上記の様な時には↑の関数が使えそうだなと思いました。
参考
https://bellcurve.jp/statistics/course/9427.html
https://bellcurve.jp/statistics/course/9936.html