仮説と推定

統計と確率の理論を何に用いるのでしょうか。データサイエンスの部分には、データとそれを生成するプロセスに関する仮説を立て、検定を行う作業が含まれます。

統計的仮説検定

データサイエンスでは、ある仮説が真である可能性が高いか否かをしめさなければならない場面がしばしば登場します。

ここで言う仮説とは「このコインには湯が編みがない」、「データサイエンティストはRよりもPythonを好む」、「鬱陶しい広告がポップアップされる上に、広告のクローズボタンが小さくて見えにくいようなWebページは内容が全く読まれることもなくページが閉じられてしまう可能性が高い」など、データに関する統計量に言い換えられる、ある種の主張のことです。統計量とは、様々な過程のもとで既知のぶんさんに従う確率変数の観測結果であると考えられ、それらの過程がどの程度正しいを与えることができます。

古典的な設定では、帰無仮説H0が基本的な立ち位置を示し、対立仮説H1と比較されます。統計量を用いて、この帰無仮説H0を棄却できるか否かをかていします。例えをしめします。

例) コイン投げ

コインに歪みがあるか否かを確認します。コインを投げて表が出る確率をpとした場合、コインに歪みがないことを締める帰無仮説は p=0.5となります。この仮設と対立仮説であるP≠0.5と比較して検定を行います。

この検定ではコインをn回投げてオモテが出た回数Xを数えます。各コイン投げはベルヌーイ思考で、XはBinomial(n, p)の確率変数となります。これは正規分布で近似可能です。

def nomal_approximation_to_binomial(n, p):
#Binomial(n, p)に相当するµとσを計算する
mu = p * n
sigma = math.sqrt(p * (1 - p) * n)
return mu, sigma

確率変数が正規分布に従う限り、実際の値が特定の範囲内に入る(もしくは入らない)確率はnomal_cdfを使って把握できます。

#変数が閾値を下回る確率はnormal_cdfで表せる
normal_probability_below = normal_cdf
#閾値を下回っていなければ、閾値より上にある
def normal_probability_above(lo, mu=0, sigma=1):
return 1 - normal_cdf(hi, mu, sigma)
#hiより小さく、loより大きければ、値はその間にある
def normal_probability_between(lo, hi, mu=0, sigma=0):
return normal_cdf(hi, mu, sigma) - normal_cdf(lo, mu, sigma)
#間になければ、範囲外にある
def normal_probability_outside(lo, hi, mu=0, sigma=1):
return 1 - normal_probability_between(lo, hi, mu, sigma)

同じことは逆の手順でも可能です。あるレベルまでの可能性の相当する区間または平均を中心とした左右対称な領域を求めます。例えば、平均を中心として60%の確率で発生する領域を求めたければ、上下それぞれの確率が20%の分を取り除けばよいです(残りは60%となります)

def normal_upper_bound(probability, mu=0, sigma=1):
#確率P(Z <= z)となるzを返す
return inverse_normal_cdf(probability, mu, sigma)
def normal_lower_bound(probability, mu=0, sigma=1):
#確率P(Z >= z)となるzを返す
return inverse_normal_cdf(1 - probability, mu, sigma)
def normal_two_sided_bounds(probability, mu=0, sigma=1):
#指定された確率を方含する(平均を中心に)対称な境界を返す
tail_probability = (1-probability) / 2
#上側の境界はている確率(tail_probability)分上に
upper_bound = normal_lower_bound(tail_probability, mu, sigma)
#下側の境界はている確率分下に
lower_bound = normal_upper_bound(tail_probability, mu, sigma)
return lower_bound, upper_bound

コイン投げの回数をn = 1000としましょう。コインに歪みはないというかせるが真であるなら、Xは平均が500で分散15.8の正規分布で近似できます。

mu_0, sigma_0 = normal_approximation_to_binomial(1000, 0.5)

実際には真であるにもかかわらずH0を棄却してしますという第一種の過誤(偽陽性)をどの程度受け入れるか、いわゆる有意性について決めておかなければなりません。年代記によると多くの場合で有意水準を5%か1%に設定しています。ここでは5%を使用します。

Xが以下d与えられる区間外になってしまったため、H0が棄却されるという状況について考えてみましょう。

normal_two_sided_bounds(0.95, mu_0, sigma=0) #(469, 531)

pが実査に0.5に等しい(つまり、H0が真である)場合を考える、Xがこの区間外となる可能性は5%であり、これは当初考えた有意性と正確に等しい値です。別の表現をしてみましょう。もしH0が偽であるにおかかわらずH0を棄却しないという第二種の過誤がおきない確率についても考えてみましょう。検定力を測るために、H0が偽と知っていることは、Xの分布に関する情報をほとんどもたらしません)コインの表が出やすいように少しだけ歪んでいて、pが実際には0.55であった場合に何が起きるのかを確認します。

この場合、検定力は次のように計算できます。

#pが0.55であると想定の元で、95%の境界を確認する
lo, hi = normal_two_sided_bounds(0.95, mu_0, sigma_0)
#p = 0.55であった場合のµとσを計算する
mu_1, sigma_1=nomal_approximation_to_binomial(1000, 0.55)
#第二種過誤とは、帰無仮説を棄却しないという誤りがあり、Xga当初想定の領域に入っている場合に生じる
type_2_probability = nomal_probability_between(lo, hi, mu_1, sigma_1)
power = 1 - typw_2_probability #0.887

代わりに帰無仮説が、コインにゆがみがない、もしくはp≤0.5であると仮定してみます。この場合片側検定を使います。Xが500よりずっと大きければ帰無仮説と棄却し、500よりも小さければ棄却しません。つまり、5%の有意性で検定を行うにはnomal_probability_belowを使って確率が95%となるカットオフ値を求めることになります。

hi = normal_upper_bound(0.95, mu_0, sigma_0) #526
type_2_probability = normal_probability_below(hi, mu_1, sigma_1)
power = 1 - typw_2_probability #0.936

Xが469より小さい(H1が真であるなら、ほとんど起こり得ない値)場合にはH0を棄却せず、一方でXが526と531の間(H1が真であるなら、多少起こり得る可能性のある値)の場合にはH0を棄却することになるため、より強い検定であるといえます。

この検定を図る別の尺度がp値である。特定の確率でカットオフを選ぶ代わりにH0が真であると仮定して、実際に観測された値を少なくとも同等に極端な値が生じる確率を計算します。コインの歪みに関する両側検定は次のように行います。

def two_sided_p_value(x, mu=0, sigma=1):
if x >= mu:
#xが平均より大きい場合、テイル確立はxより大きい分
return 2 * normal_probability_above(x, mu, sigma)
else:
# xが平均より小さい場合、テイル確率はxより小さい分
return 2 * normal_probability_below(x, mu, sigma)

オモテが530回出た場合は次のように計算できます。

two_sided_p_value(529, mu_0, sigma_0) #0.062

これが理にかなった推定であることを納得するために、シュミレーションを行います。

extreme_value_count = 0
for _ in range(100000):
num_heads = sum(1 if random.random() < 0.5 else 0
#1000回のコイン投げを行いオモテが出る回数を数える
for _ in range(1000))
if num_heads >= 530 or num_heads <= 470:
extreme_value_count += 1
#そのうち極端な回数はどれだけ出たかを数える
print(extreme_value_count / 100000) #0.062

5%よりも小さい値となたので、帰無仮説を棄却します。これは先に行った検定とつぎも同じです。

upper_p_value = normal_probability_above
lower_p_value = normal_probability_below

片側検定の場合、525回のオモテが出れば、帰無仮説を棄却しませんが、

upper_p_value(524.5, mu_0, sigma_0) #0.061

527回ならば、棄却することになります。

upper_p_value(526.5, mu_0, sigma_0) #0.047

信頼区間

分布を未知のパラメータとして、コインのオモテが出る確率に関する仮説検定を行ってきました。もしこれが本当であるなら、観測値の周辺の信頼区間をもとめるのが3番目の手法となります。

例えば、オモテを1、ウラを0とするベルヌーイ変数の平均を見ることで歪みのあるコインに対する確率を推定できます。1000回の試行で525回のオモテが出たとすると、pの推定値は0.525です。

この推定値はどの程度信頼できるでしょうか。pの正確な値を知っているなら中心極限定理により、このベルヌーイ変数の平均値は近似的に平均p及び次の標準偏差の正規分布に従います。

math.sqrt(p * (1 - p) / 1000

ここではpが未知となっているので、先程の推定値をつかいます。

p_hat = 525 / 1000
mu = p_hat
sigma = math.sqrt(p_hat * (1 - p_hat) / 1000) #0.0158

この値は完全に理にかなったものだとはいえませんが、ともかくこの方法が使われます。正規分布の近似をつかうと、pの正しい値が次の区間に入るのは「95%の確率で信頼できる」という結論になります。

normal_two_sided_bounds(0.95, mu, sigma)#[0.4940, 0.5560]

この結果、0.5はこの信頼区間内にあるため、コインに歪みがあるとは結論付けられません。では、重手が540回出た場合はどうでしょうか。

p_hat = 540 / 1000
mu = p_hat
sigma = math.sqrt(p_hat * (1 - p_hat) / 1000) #0.0158
normal_two_sided_bounds(0.95, mu, sigma) #[0.5091, 0.5709]

コインに歪みがないとした場合の値は、この信頼区間内にはまってません(つまり仮説が正しいとすれば95%の確率でその範囲に入るという検定に対して、この「コインに歪みがない」という仮説は成立しません)

pハッキング

5%の確率で誤って帰無仮説を棄却する手順は、定義により5%の確率で誤って帰無仮説を棄却します。

def run_experiment():
#歪みのないコインを1000回なげて、オモテが出たらTrue,ウラが出たらFalse
return [random.random() < 0.5 for _ in range(1000)]
def reject_fairness(experiment):
#5%の有意水準を用いる
num_heads = len([flip for flip in experiment if flip])
return num_heads < 469 or num_heads > 531
random.seed(0)
experiments = [run_experimet() for _ in range(1000)]
num_rejections = len([experiment in experiments if reject_fairness(experiment)])
print num_rejections #46

これが意味するのは、「有意な」結果を得ようとすれば、それは可能だということです。エータセットに対する十分な仮説検定を行えば、そのうちの1つは明らかな有意性を示します。外れ値を適切に取り除くことで、p値はおそらく0.05未満にできます。(相関でも同じことができます)

「p値を使った推定の枠組み」から得られる結論に対して何らかの手を入れてしまうこの手法をpハッキングと呼びます。この手法を批判した優れたThe earth is Roundという記事があります。

適切なデータサイエンスを行いたいのであれば、データを調査する前に仮説を決定し、データの整理は仮説を前提とせずに行い、p値は常識の代用品とはならないことに注意しないといけません(これと異なる手法をベイズ推定で、後ほど扱います)

事例:A/Bテストの実施

ユーザ体験の最適化、直接的に言うとユーザに広告をクリックさせるのがデータサイエンティストにとって重要な任務の1つとなります。

例えば、データサイエンティスト向けの新しい栄養ドリンクを開発しました。そこで、広告A「最高に美味しい」と広告B「こんなにバイアスが減った」のどちらがいいかを決めます。

Aを見た1000人のうち990人が広告をクリックしたのに対し、Bは10人しかクリックしなかったとすれば、Aの方がよいとわかりますが、違いがそれほど明確でなければどうでしょうか。そこで統計的推定が役立ちます。

Na人が広告Aをみて、そのうちのna人がクリックしたとします。これは、拘束が表示される確率がpaのベルヌーイ試行とみなせます。(Naが十分におおきいとして)na/Naは、平均pa、標準偏差σa=(pa(1-pa)/Na)**0.5の正規分布で近似できる確率変数となります。
同様に、nb/Nbは、平均pbで標準偏差σb = (pb(1 — pb) / Nb)**0.5の正規分布で珍事で確率変数となります。

def estimated_parameters(N, n):
p = n / N
sigma = math.sqrt(p * (1 - p) / N)
return p, sigma

これらの正規分布が独立であるなら(それぞれの別のベルヌーイ試行であるはずなので、この前提は妥当だと思われます)、その差も正規分布となり、平均pb-paおよび標準偏差(σa**2 + σb**2)**0.5となるはずです。

これにより、標準正規分布を持つ統計量を用いて、paとpbが等しいという帰無仮説を検定できます

def a_b_test_statistic(N_A, n_A, N_B, n_B):
p_A, sigma_A = estimated_parameters(N_A, n_A)
p_B, sigma_B = estimated_parameters(N_B, n_B)
return (p_B - p_A) / math.aqrt(sigma_A ** 2 + sigma_B ** 2)

例えば「最高においしい」をクリックしたのが1000人中200人。「こんなにバイアルが減った」をクリックしたのが1000人中180人だった場合、次のように計算できます。

z = a_b_test_statistic(1000, 200, 1000, 180) #-1.14

平均が等しいときに、この大きさの違いが生じる確率は次の値となります。

two_sided_p_value(z) #0.254

これは値として十分に大きいため、帰無仮説を棄却できず、違いがあるという結論は導けません。一方「こんなにバイアスが減った」をクリックしたのが1000人中150人だった場合、

z = a_b_test_statistic(1000, 200, 1000, 150) # -2.94
two_sided_p_value(z) # 0.003

両方の広告が等しい場合に、このクリック数の違いが出る確率は0.003しかないということになります。

ベイズ推定

先の手法は、検定の結果「帰無仮説が正しい場合に、これだけの極端な違いが出る確率は0.3%しかない」という内容を導いています。
これとは別に、未知のパラメータを確率変数として扱う推定の手法があります。アナリストはパラメータの事前分布に対して観測データとベイズの定理を用いて、パラメータの事後分布を求めます。検定結果の確率を使う代わりに、パラメータ自身を用いて確率的な判断を行います。

例えば、(コイン投げなど)確率が未知のパラメータであったとすると、0か1の間の様々な値を取りうる分布が事前分布として頻繁に使われます。

def B(alpha, beta):
#確率の総和が1となるように定数で正規化する
return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)
def beta_pdf(x, alpha, beta):
if x<0 or x>1: #[0,1]の間では重みは0となる
return 0
return x ** (alpha - 1) * (1 - x) ** (beta - 1)/B(alpha, beta)

一般的には、この分布は次の重みを中心とします。

alpha / (alpha + beta)

そして、alpha,betaが大きければ、分散は狭くなります。
例えば、alpha,betaが共に1であった場合、単なる一様分布となります(中心の0.5から一様に分散します)betaよりもずっとalphaが大きければ、1の布巾に重みが偏ります。逆にbetaよりもずっとalphaが小さければ、重みは0の知覚に位置します。

pの事前分析について考えましょう。コインに歪みがある否かを明言したくないなら、alpha,betaを共に1とします。または、55%の確率で表が出るという強い確信があるなら、alpha=55, beta=45にします。

続いてコイン投げを何度も行い、表が出た回数と裏が出た回数tを観察します。ベイズの定理(くわえて、この話題wお読み進むための多少退屈な数学)によると、pの事後分布もβ分布となりますが、そのパラメータは、alpha + h, beta + tとなります。事後分布もβ分布となれるのは、偶然の一致ではなく、オモテの出る確率は二項分布により与えられ、ベータ分布は二項分布の共役事前分布です。つまり、二項分布の観測データによりベータ分布の事前分布を更新するとベータ分布の事後分布が得られます。

ここで、10回コインを投げて3回しかオモテが出なかった場合、一様の事前分布(ある意味、コインに歪みが歩かないか、明確にしていない)から始めたとしても、事後分布は0.33を中心としたBeta(4, 8)になるでしょう。はじめは度の確率も等しく起こりうると考えていたので、推測値は観測された確率に近くなります。

Beta(20, 20)で(コインにはおよそ歪みがないと表明して)始めた場合、事後分布は0.46を中心としたBeta(23, 27)になるでしょう。これは裏が多く出るようなバイアスを持っていると更新されたことを示します。

事前分布をBeta(30,10)とした場合、(75%の割合でオモテが出ると考えている)事後分布は0.66を中心にBeta(33, 17)となるでしょう。この場合、引き続きオモテが出るバイアスをもっていると考えられますが、当初の予想よりはその偏りは小さい値に更新されたと思います。

コイン投げの試行回数をもっと増やせば、事前分布の影響力は低下し、最終的には不全分布がどうであっても同じ事後分布に近づくでしょう。

例えば、事前にコインのバイアスをどのように考えていたとしても、コインを2000回投げた内オモテが1000回出るのを見れば、その考えを維持するのはこんなんです。ここで興味深いのは仮説の確率に関して「事前分布と観察されたデータからコインの表が出る確率が49%と51%の間に存在する可能性は、5%しかない」といえることです。これは「コインに歪みがないとして、極端な値がでる可能性は5%しかない」と言及するのとは大きな違いがあります。

ベイズ推定を余地いて行う仮説検定は、しばしば議論の対象になります。使用する数学が多少複雑であることが一因ですし、事前分布の選択が主観的である点も一因です。

Like what you read? Give Okazawa Ryusuke a round of applause.

From a quick cheer to a standing ovation, clap to show how much you enjoyed this story.