はじめまして、ZOZO研究所福岡の家富です。画像検索システムのインフラ、機械学習まわりを担当しています。

今回は、t検定におけるサンプルサイズが与える影響を解説します。

t検定の使われ方

近年、施策が有効かどうかをデータを元に統計的に判断していこう、という話を聞くことが増えてきました。

経済学の流行においても、統計的な指標を重要視する流れが強まってきています。例えば、貧困対策にお金をどの程度どのような用途で支給するのが良いか、といった議論で利用されることも多くなってきています。

www.amazon.co.jp

Web業界においても、サイトの変更や施策の有効性をA/Bテストなどを実施し、統計的に判断していく流れが主流になってきています。なお、本記事では以下の文献を「A/Bテスト本」と呼びます。

www.amazon.co.jp

そして、このようなA/Bテストから得られた統計情報に対し、それが有効か否かを判断する方法の1つにt検定があります。

t検定

例として、次のような各ユーザーの一週間の購入金額のデータを考えます。

2つの施策AとBがあり、それぞれの施策に対し、以下のような購入者の購入金額のデータが得られたとします。

施策Aの購入者のデータ： $\{ X_1, X_2,...,X_n \}$
同様に施策Bの購入者のデータ： $\{ Y_1, Y_2,...,Y_m \}$

ここで、n、mはそれぞれ施策A、Bの購入者数を表し、 $X_i$ 、 $Y_i$ はそれぞれのユーザーの購入金額とします。

この時に、「施策Aの分布と施策Bの分布を比べて、施策Aの方が良いと言えるのか？」というのが一般的に考えたい問題かと思います。

一番最初に思いつく方法は、それぞれの平均値 $\frac{1}{n}\sum_iX_i$ 、 $\frac{1}{m}\sum_iY_i$ を見る方法です。 $\frac{1}{n}\sum_iX_i > \frac{1}{m}\sum_iY_i$ ならば施策Aの方が良いと考えます。

ここで問題になるのは、「たまたま施策Aの平均の方が少し良かっただけなのでは？」という可能性を捨てきれない点です。非常に小さな差しかなければ、「たまたまじゃないのか」という気になります。では、「どれほどの差があれば十分施策Aの方が良いと言えるのか？」というの疑問が湧いてきます。検定を使えば、この疑問に答えることができます。

t検定の場合、次のような論法をとります。

まず、施策A・Bのどちらのデータも、ある潜在的な確率分布から独立にサンプリングされたものと考えます。なお、この見方は実際のデータに対してよく仮定されるものです。次に「施策Aの潜在的な分布の平均値」と「施策Bの潜在的な分布の平均値」は等しいと仮定します。そして、以下の式を計算します。

$t = \frac{ \bar{X} - \bar{Y} } { \sqrt{ \frac{s(X)^{2}}{n} + \frac{s(Y)^{2}}{m} } }$

この、 $\bar{X}, s(X)^{2}$ は以下で定義されるものです。 $\bar{X} = \frac{1}{n}\sum_iX_i, s(X)^{2} = \frac{1}{n}\sum_i(X_i - \bar{X})^{2}$

「ある仮定のもと（詳細は後述します）」で、この値は自由度m+n–2のstudentのt分布に従うことが証明されています。正確には、ウェルチのt検定なので自由度が多少異なりますが、m=nの時はほぼ同様の自由度となります。

en.wikipedia.org

そのため、t分布の表を見ることで、上記で得られたt値が実際にどの程度の確率で生じたものなのかがわかります。

例として、m=n=121の場合、つまり該当のt分布の自由度が240となる場合を考えてみます。計算されたt値が1.392であったとします。これは自由度240のt分布において上位5％となる点は1.651なので、そこには含まれません。一方、上位10％でのt分布の値は1.285です。そのため、上記の値は5〜10％の間で生じる値ということになります。

検定でよく採用される5％を基準とした場合、上記の値は起こりうる確率の範囲と捉えることができ、特に矛盾はないと考えます。そして「施策Aの潜在的な分布の平均値」と「施策Bの潜在的な分布の平均値」が等しいことに問題はないと解釈します。

また、10％を基準とした場合には、上記の値は起こりうる確率の範囲外にあると考えられます。すなわち、「施策Aの潜在的な分布の平均値」と「施策Bの潜在的な分布の平均値」は等しくない、施策Aは有意に有効であると考えられます。

上述のように、2つのデータの違いを評価することに使われます。

この「破られることを期待される仮定を置く」という論法が若干わかりにくくしている部分であり、時々誤解を生んでいることもあります。このような論法が使われる背景には、たくさん施策し数値的に特に目立った施策だけを拾っていこう、という発想があります。

t検定の問題点

前章でt検定の使用方法を説明しました。実際に使用する際には、t値がt分布に従うかどうかがポイントになります。

「ある仮定のもと」と前章で述べましたが、ここではその点をさらに見ていきます。t値がt分布に収束するには、以下のように定義された値 $t_{prev}$ がt分布に従う必要があります。

$t_{prev} = \frac{\bar{X} - \mu}{ \sqrt{ \frac{ s(X)^{2}}{n}} }$

なお、 $\mu$ は「潜在的な分布の平均値パラメータ」であり、データから推測はできますが、一般のデータからは得られないことに注意が必要です。

次に、 $t_{prev}$ がt分布に従うための条件を見るために、以下のように分解します。

$t_{prev} = \frac{\sqrt{n}(\bar{X} - \mu)}{\sqrt{s(X)^{2}}} = \frac{\frac{\sqrt{n}(\bar{X} - \mu)}{\sigma}}{\sqrt{\frac{s(X)^{2}}{\sigma^{2}}}}$

なお、 $\sigma$ は「潜在的な分布の標準偏差パラメータ」であり、データから推測はできますが、一般のデータからは得られないことに注意が必要です。

ここでは、 $\frac{\sqrt{n}(\bar{X} - \mu)}{\sigma}$ が中心極限定理より「nが十分大きい時に」正規分布に収束します。また、 $\frac{ns(X)^{2}}{\sigma^{2}}$ が「nが十分大きい時に」自由度n-1のカイ二乗分布に収束します。このことから、 $t_{prev}$ は「nが十分大きい時に」自由度n-1のt分布に従うと言えます。なお、nが十分大きい時なので、nとn-1の差はほぼ考えなくても良いです。

この時「nが十分大きい時に」という表現が、実際にどれくらいの大きさならば誤差を無視できるのかが問題となります。

「A/Bテスト本」では、元の分布の歪度（ $\beta_1(X)$ ）が大きい場合の問題を指摘しています。なお、 $\beta_1(X)$ は以下のように定義されます。

$\beta_1(X) = \frac{E( X - E(X))^{3}}{Var(X)^{3/2}}$

なお、 $E(X)$ は分布の期待値、 $Var(X) = E( X - E(X))^{2}$ で定義される分布の分散を表します。歪度が大きい場合、中心極限定理による近似度があまり良くありません。そのため、t検定をする際にはサンプルサイズを大きく取る必要があると主張しています。

「A/Bテスト本」では、サンプルサイズnは $355 \times \beta_1(X)^{2}$ 以上必要だと主張しており、歪度が10以上となるとかなりのサンプルサイズが必要だと主張しています。さらに、この主張がどこからくるのかについて調べていくと「How Large Does n Have to Be for Z and t Intervals?, Dennis D. Boos and Jacqueline M. Hughes-Oliver」という論文によるものだとわかりました。

次章では、上記論文の内容を紹介していきます。

論文手法

「How Large Does n Have to Be for Z and t Intervals?, Dennis D. Boos and Jacqueline M. Hughes-Oliver」（以下「本論文」と呼ぶ）ではまず、歪度と中心極限定理の関係について述べています。

1次元の実数上の確率分布に対し、正規分布との関係を述べた定理としてGram–Charlier A seriesがあります。この定理を $\frac{\sqrt{n}(\bar{X} - \mu)}{\sigma}$ に対して適用したものがEdgeworth seriesです。en.wikipedia.org

この定理により、一般の確率分布の累積分布関数をn^-1/2の級数展開として得ることができます。そして、n^-1/2以降の項はn^-1,n^-3/2と続きます。サンプルサイズnが大きくなる場合、n^-1以降は収束が速いため、n^-1/2の項のみに注目すれば良いことがわかります。ここで、n^-1/2の項の係数をみると歪度が出てくるため、収束の速度に歪度が重要なパラメータとして関わっていることもわかります。また、中心極限定理の収束速度を表す定理の1つであるBerry-Esseenの定理においても、やはり収束誤差を $n^{-1/2} \times \beta_1(X)$ で抑えられると言われています。

en.wikipedia.org

上記の知見を元にし、本論文では「上位5%を表す閾値0.05でのt値に対し、 $c \times n^{-1/2} \times \beta_1(X)$ だけずれる。なお、cは分布に依らない定数。」と仮定します。一般にこの閾値は「 $\alpha=0.05$ 」のように $\alpha$ で表すことが多いため、本記事でもこの表記を使用します。

次に、先程登場したcを求める方法について述べます。本論文では様々な既存の分布から実際にサンプリングを行い、歪度とズレに関して線形回帰をすることによって求めています。よく知られた解析的な確率分布を用いているので、 $\mu$ のような潜在的なパラメータをシミュレーション結果からではなく、実際に計算できるので、結果としてt値を計算することが可能となります。

本論文では以下の分布を使用しています。

ガンマ分布
- $\theta$ = 1, k=1, 1.78, 4, 16としたもの

ja.wikipedia.org

ワイブル分布
- scale parameter=1, shape parameter= 1.2, 1.6, 2.2としたもの

en.wikipedia.org

論文「An Asymptotically Distribution-Free Test for Symmetry versus Asymmetry」による3パターンのオリジナル分布
- lam1, lam2, lam3, lam4 = 0, 1.0, 1.4, 0.25 # 7
- lam1, lam2, lam3, lam4 = 0, 1.0, 0.000070, 0.10 # 8
- lam1, lam2, lam3, lam4 = 0, -1.0, -0.1, -0.18 # 12

上記の分布に対し、本論文においてはズレ（miss）を次のように計算します。

まず、各分布に対して30サンプル抽出した時のt値を10,000セット用意します。なお、今回の場合だと300,000サンプル合計抽出します。

次に、用意した10,000個のt値を大きい順に並べ替え、大きい方からt値を調べていきます。正確な自由度30のt分布においては5％に相当するt値=1.697と比較していき、それを下回るのが何番目になるかをチェックします。例えば、上位から125番目が該当した場合、上位5％とのずれは0.05 - 125/10000 = 0.0375となります。

このように様々な歪度の分布からのt値を計算し、実際のズレのデータを取得します。論文では、さらに小さい方からのt値も見ていき、左側検定の場合のズレも計算しています。

以下はズレを取得する処理のPythonの疑似コードです。

t_list = []
for s_id in range(10000):
    np.random.seed(s_id)
    xs = []    
    for _ in range(30):
        # scale parameter=1, shape parameter=1 のgamma分布からのサンプリング
        a = np.random.gamma(1, 1)
        xs.append(a)
    t_list.append( get_t(xs) )  # t値を計算して、リストに登録
    
t_list.sort()
for i,t in enumerate(t_list): # left
    if t > - 1.697:
        left_percentile = i/10000
        break
for i,t in enumerate(reversed(t_list)): # right
    if t < 1.697:
        right_percentile = i/10000
        break
print( 0.05 - left_percentile, 0.05 - right_percentile)

本論文では、左側のズレが0.116、右側のズレが-0.067という結果が掲載されています。一方で、私が行った実験結果では左側のズレが0.11186489、右側のズレが-0.03688539という結果が出ています。

最後に、この結果を用いて必要なサンプルサイズを求めます。本論文ではズレ（miss）が $\alpha$ =0.01以下となるようにサンプルサイズnを決めます。

$miss = c \times n^{-1/2} \times \beta_1(X)$

これをnの式に直すと、 $n = (\frac{c}{miss})^{2} \times \beta_1(X)^{2}$ となります。

先程の実験結果から、値の大きい左側を使うと以下の結果が算出されます。

本論文結果： $134.56 \times \beta_1(X)^{2}$
私の実験結果： $125.13 \times \beta_1(X)^{2}$

誤差を0.005以下にすると考えた場合、係数は本論文では538.24、私の実験結果では500.55となります。「A/Bテスト本」で述べられていた355という係数は得られませんでしたが、数値のオーダーはほぼ同等のものが得られたと考えられます。

実際の購入金額データに対する考察

冒頭で紹介した「各ユーザーの一週間の購入金額のデータ」に対して歪度を計算してみました。歪度は700.43となり、精度を得るために必要なサンプルサイズはn=61,389,051.39と算出されます。これは現実的ではない値です。

「A/Bテスト本」ではこのような場合、ある閾値以上の部分は落とすという方法を推奨しています。今回は10万円以上のユーザーは、別データに分けました。そうすることで、10万円未満のデータの分布は歪度が2.750となり、現実的な数値となりました。

なお、10万円以上のデータの分布はこれにより歪度が35.637となり、だいぶ小さい値になりました。扱いやすくなりましたが、まだ工夫が必要な値です。さらに細かく層を分けて対応するなどの工夫が必要でしょう。

まとめ

本記事では、t検定におけるサンプルサイズの考察を紹介しました。

実際にt検定の理論を深堀りした際には、Edgeworth seriesなどは意外と文献が見つけにくく、苦労しました。しかし、自ら数値実験をすることで実際の数値とのズレ具合など、知見を得ることができました。また、実データは想像以上に歪みの大きいデータが多いため、そのまま適用とはいかないケースがかなりあることも実感しました。