チョコボール統計

チョコボールの秘密を統計解析で明らかにしていく。おもちゃのカンヅメ欲しい。

エンゼルを高確率で当てられる運の良さを評価(階層ベイズモデルによる個人差の推定)

【概要】

  • 恒常的に「運」が良いなんて本当にあるのでしょうか?
  • 銀のエンゼルの出現確率が、チョコボールの購入者によって差異が出るのかを調査することで、運の良さを定量的に評価してみます。
  • 購入者毎に独立に出現確率を推定、一般化線形混合モデルを利用して”個人差”として「運」を推定、階層ベイズモデルで購入者毎の出現確率を推定と三つのモデルで評価してみます。
  • 結果は、今回のデータでは「運」の大きさに統計的な差異があるとは言えないということがわかりました。

f:id:hippy-hikky:20190309151311j:plain

【目次】

はじめに

「あの人は運が良いよね」とか「俺は運が悪い」っていうセリフ、聞いたり言ったりすることよくありますよね?
でも、恒常的に「運」が良いなんて本当にあるのでしょうか? 本当はランダムな現象であるにも関わらず、認知バイアスにより成功(失敗)の記憶が強化されているだけではないでしょうか??

ということで今回は、 チョコボール銀のエンゼルの出やすさ(出現確率)が購入者によって差があるのかを分析し、 購入者の運の強さを測れるのかを試してみます。

もし統計的に明らかに発生確率に差が生じているのなら、 それは本当に神秘的な何かか、不正が行われているのかのどちらかのはず。

(内容に間違いや気になった点がありましたら、ご指摘くださると助かります)

【トップに戻る】

問題設定

当ブログでは、チョコボールを買ってきて計測することをライフワークとしています。 これまでにもちょくちょく言及していますが、 チョコボールの購入者は筆者のみではありません。 友人、同僚から寄付(チョコボール)をいただくことがよくあります。

その中でも、よく寄付してくれる同僚が二人います(同僚B,Cと呼ぶ、本記事の中でAは筆者を指す)。 同僚Bは人格者で周りから好かれている良いやつです。 一方、同僚Cはちょっとクセがあります。 この二人、どちらも頻繁に寄付をしてくれるのですが、 同僚Bからもらったチョコボールから銀のエンゼルがやたらと出る感じがします。 (一方Cは。。。以下略)

そこで、筆者(A)、同僚B、同僚Cが購入したチョコボールのエンゼル出現確率に差があるのかを調査してみます。 ランダムにチョコボールを購入しているのであれば、統計的に差異は生じないはずです。 本当に差があるのだとしたら、それは神秘的な作用によるものか、 または、不正(エンゼル出現に関するなんらかの情報を持っている)が行われているのか。。。

【トップに戻る】

アプローチ

チョコボールのエンゼル出現確率の推定は、ベイズ推定を利用します。 エンゼル出現率とベイズ推定については当ブログの以下の記事を参照ください。*1

ベイズ推定でエンゼルの出現確率を予測する - チョコボール統計

ベイズ推定で【銀の】エンゼルの出現確率を予測する - チョコボール統計

チョコボールに学ぶ実践的ベイズ統計モデリング入門 #1 - チョコボール統計

今回は、上記の記事のモデルをそのまま拡張して、 購入者毎に独立に出現確率を推定してみます。 その後モデルを拡張し、個人差を加えた階層モデルとしてエンゼルの出現確率を推定し、 比較を行ってみます。 それぞれのモデルの詳細は後述します。

【トップに戻る】

データ

利用するデータは当ブログで集計しているデータです。 全データから、「金のエンゼル2倍キャンペーン」の期間のデータを除外します。 このキャンペーンでは銀のエンゼルが含まれないと公式で発表されているからです。

データの概要

  • データの集計期間:2017/11/19〜2019/03/07
購入者 購入数 銀のエンゼル 頻度
A 142 6 0.042
B 80 5 0.063
C 59 2 0.034
total 281 13 0.046

全体での最尤推定量(頻度)は4.6%となっています。

【トップに戻る】

購入者毎のエンゼル出現確率の推定

チョコボールはランダムに購入すると仮定し、銀のエンゼルが出るか出ないかはある一定の確率に従って 毎回ランダムに決まるとします(つまりベルヌーイ試行)。 そのため、基本的なモデルとして銀のエンゼルが出る確率を\thetaとした二項分布を考えます。

{ \displaystyle
p(y | \theta) = \binom{n}{y}\theta^{y}(1-\theta)^{n-y}
}

エンゼル出現確率の推定に利用したコードはひとまとめにして5.5に貼り付けておきます。 参考にしていただけたらと思います。

最尤推定

まずは最尤推定量を確認しましょう。 最尤推定については、以下の記事などを参照ください。

チョコボールに学ぶ実践的ベイズ統計モデリング入門 #1 - チョコボール統計

データの性質を理解する - 機械と学習する

上記の記事などを参照すると、二項分布の最尤推定量は頻度であることがわかるので、 4.1節の表に記載の通りです(図で示すと以下の通り)。

f:id:hippy-hikky:20190309154627p:plain
購入者A,B,C毎のエンゼル出現確率の最尤推定

これだけ見ると、Bの運が良いように見えると思います。 しかし、この差は本当に差があるのでしょうか?偶然とは考えられないでしょうか?

これは点推定なので、この差がどのくらい意味があるかはまだわかりません。 検定などすれば差の評価はできますが、 以降の節ではベイズ推定を使ってエンゼルが出る確率\thetaの 事後分布を利用して評価していきます。

ソースコード解説(クリックで展開) 最尤推定量の計算とプロットは、5.5のセル番号6です。
最尤推定量は頻度なので、単純にエンゼルの数を全体数で割っているだけなので、解説は不要かと思います。

独立に出現確率を推定

ベイズ推定で\theta_i(iは購入者を指すindex)を推定しましょう。

モデルは、購入者毎の購入数n_iとエンゼル出現確率\theta_iをパラメータに持った二項分布です。

{ \displaystyle
p(y_i | \theta_i) = \binom{n_i}{y_i}\theta_i^{y_i}(1-\theta_i)^{n_i-y_i}
}

ベイズ推定なので事前分布を仮定する必要があります。 事前分布にはベータ分布(下記)を採用します。

{ \displaystyle
p(\theta) = Beta(\alpha, \beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}
}

ベータ分布は二項分布の共役事前分布ですし、 \alpha\betaの二つのパラメータで様々な形状を表現できます。

このモデルを図で表すと以下のようになります(手書きでごめんなさい)。

f:id:hippy-hikky:20190309160816p:plain:w300
購入者毎に独立にエンゼル出現確率推定のモデル

\thetaの事前分布であるベータ分布のパラメータを(\alpha, \beta)=(1,1)、 つまり無情報事前分布として推定すると、結果は以下のようになりました。

f:id:hippy-hikky:20190309161227p:plain
購入者毎に独立に出現確率を推定した結果(\thetaの事後分布)

この結果から、MAP(事後確率最大の点)は5.1に記載の最尤推定とほぼ同じように見えますが、 それぞれの分布は重なっており、購入者毎の明確な差は見えなそうだといえます。

事後分布の差の分布も見てみましょう。

f:id:hippy-hikky:20190309161708p:plain
事後分布の差の分布。左からA対B、A対C、B対C。

この図は、購入者二人ずつの事後分布の差の分布です。 0付近に分布が集中していれば、差はないと言えます。 この図から、B対Cで少し差の分布が0から離れていますが、95%信用区間の中に収まっていますので、 差があるとは言えないということがわかります。

ソースコード解説(クリックで展開) 購入者毎に独立にエンゼルの出現確率を推定している部分は、5.5のセル番号7です。

with pm.Model() as comparing_buyer_i:
    theta = pm.Beta('theta', alpha=1, beta=1, shape=len(set(buyer_idx)))
    angel = pm.Binomial('angel', n=total_counts[lst_buyer], p=theta[lst_buyer], observed=angel_counts[lst_buyer])
    trace_iso = pm.sample(5000, chains=1, random_seed=100)

pymc3の記法として、withブロックで囲んだ中でモデルを定義します。 上から、\thetaの事前分布として無情報事前分布を設定し、 二行目で尤度として二項分布を用いることにしています。 三行目で、乱数のシードを固定してMCMCサンプル5000取得しています。

8~10セルまでは、MCMCサンプルを可視化して事後分布を確認しています。 この中では、10セル目で、購入者毎の差の分布を可視化していますが、 事後分布をMCMCサンプルで得ているので、単純に差をとった分布を確認すれば良いというのが お手軽で良いですね。

購入者の個人差をパラメータに加えて推定

5.2では、購入者独立にエンゼル出現確率\theta_iを推定しました。 5.2の結果だけでも、統計的に明確な差は見えないということがわかりましたが、 もう少しモデルを複雑化してみましょう*2

エンゼルの出現確率は、森永さんの設定により、ある決まった値になっているものと仮定します(これをベースの出現確率と呼ぶ)。 もし人によって「運」の量のような物があるのなら、 ベースの出現確率に個人毎が持っている運の要素が加えられて調整されていると考えます。

これをモデル化してみましょう。

f:id:hippy-hikky:20190309163347p:plain:w300
購入者の個人差(運)を加えたモデル

\thetaは確率なので、[0, 1]の範囲に収まる必要があります。 そのため、ロジット関数を使ってエンゼル出現の対数オッズがp+u_iで表されると仮定します。 ここで、pは森永さんが設定するベースの出現確率を決めるパラメータで、 u_iは購入者iの運の強さを表します。

pu_iに適当な事前分布(裾が広い一様分布状の正規分布)を仮定して推定した結果は以下の通りです。

f:id:hippy-hikky:20190309164118p:plain
購入者毎の運を定めるパラメータu_iの事後分布

この結果から、運の要素というのは非常に小さいということがわかりますし、購入者毎の差もほとんどなさそうなことがわかります。 u_iの事後分布の差の分布を以下に示します。

f:id:hippy-hikky:20190309164357p:plain
運の要素の事後分布の差の分布。左からA対B、A対C、B対C。

B対Cを見ると多少差が大きいようですが、信用区間の中に収まっており、明確に差があるとは言えないようです。

森永さんの設定する真の出現確率については、以下のように予測できました。

f:id:hippy-hikky:20190309165721p:plain
真のエンゼル出現確率(対数オッズ)の事後分布

この図の横軸は、対数オッズ\log{\frac{q}{1-q}}です。 確率に変換するには、シグモイド関数を通す必要があります。 結果として、95%信用区間は以下の表の通りとわかりました。

2.5%点 平均 97.5%点
0.25% 4.19% 34.00%

まだ範囲が広いですね。データをもっと集める必要があります。

ソースコード解説(クリックで展開) 購入者毎の個人差を加えたモデルを計算している部分は、5.5のセル番号11です。

with pm.Model() as comparing_buyer_m1:
    su = pm.HalfNormal('su', sd=10)
    p = pm.Normal('p', mu=0, sd=10)
    u = pm.Normal('u', mu=0, sd=su, shape=len(set(buyer_idx)))
    
    angel = pm.Binomial('angel', 
                        n=total_counts[lst_buyer], 
                        p=pm.math.sigmoid(p+u[lst_buyer]), observed=angel_counts[lst_buyer])
    
    trace_h1 = pm.sample(3000, chains=3, random_seed=100)

ロジット関数を展開すると、二項分布のパラメータ\thetaシグモイド関数で表現できるので、 二項分布にはp=pm.math.sigmoid(p+u)という形で入力しています。 推定対象のパラメータはp, u_i, suの3種類ですので、 それぞれの事前分布を設定しています。 なお、suは個人差を決めるパラメータu_iの事前分布のパラメータです。

セル番号16,17では、対数オッズで予測されるパラメータを確率に変換するためにシグモイド関数を定義して、 ベースの出現確率を算出しています。

階層モデルとして予測

もう少し見方を変えて、別のモデルを立ててみましょう。 5.3節では、購入者毎のエンゼルの出現確率は、 ベースの出現確率pに個人差u_iが加わったものとしてモデルを作りました。

ここでは、森永さんの設定するベースの出現確率の確率分布から購入者毎の出現確率がサンプルされると考えます。 つまり、5.2節のように購入者毎に独立に出現確率が決まるのではなく、 ベースの確率になんらかの調整が加わって購入者毎の確率が決まると考えます*3

これは、階層ベイズモデルを考えることでモデル化できます。

モデルを以下の図に示します。

f:id:hippy-hikky:20190309173044p:plain:w250
購入者毎のエンゼル出現確率は、真の確率(事前分布)からサンプルされると仮定するモデル

階層事前分布Beta(\alpha, \beta)のそれぞれの事前分布は裾の広い正規分布(無情報事前分布)として、 上記のモデルに従って購入者毎のエンゼル出現確率を推定した結果を以下に示します。

f:id:hippy-hikky:20190309173720p:plain
階層モデルで推定した購入者毎のエンゼル出現確率の事後分布

購入者毎の差の分布は次の通りです。

f:id:hippy-hikky:20190309173848p:plain
階層モデルで推定した購入者毎のエンゼル出現確率の事後分布の差の分布

これらの図から、やはり購入者によって差はほとんど無いということがわかりました。

では、\thetaの事前分布、つまり、ベースの出現確率の分布を確認してみます。

f:id:hippy-hikky:20190310220856p:plain
\thetaの事前分布。濃い青線が期待値を示す。オレンジはMCMCサンプルからのランダムなサンプル。

この結果から、ベースの出現確率の期待値は4%ということがわかりました。

ソースコード解説(クリックで展開) 購入者毎に独立ではなく、階層ベイズモデルでエンゼルの出現確率を推定している部分は、 5.5のセル番号18です。

with pm.Model() as comparing_buyer_m2:
    alpha = pm.HalfCauchy('alpha', beta=10)
    beta = pm.HalfCauchy('beta', beta=10)
    
    theta = pm.Beta('theta', alpha=alpha, beta=beta, shape=len(set(buyer_idx)))
    
    angel = pm.Binomial('angel', n=total_counts[lst_buyer], p=theta[lst_buyer], observed=angel_counts[lst_buyer])
    
    trace_h2 = pm.sample(2000, chains=3, random_seed=100)

5.2の購入者毎に独立に出現確率を推定したコードと似ています。 5.2では、\thetaの事前分布はBeta(\alpha=1,\beta=1)と固定されていましたが、 ここでは\alpha, \betaも推定対象です。 そのため、これらのパラメータの事前分布を定義しています。

セル番号21では、\theta_iの事前分布をプロットするコードを書いています。 これは、参考文献2からほぼそのままコピーしてきました。

ソースコード

まとめ

以上、購入者毎のエンゼル出現確率に差があるのかを3つのモデルを使って推定しました。

最尤推定(5.1)だけでは、意味のあるほどの差なのかわからなかったのですが、 ベイズ推定により事後分布を推定(5.1, 5.2, 5.3)したところ、 購入者毎にエンゼルの出現確率が明らかに異なるとは言えないということがわかりました。

また、購入者毎の効果(運)を明確に分離するモデル(5.3と5.4)を使って森永さんが設定しているであろう 銀のエンゼルのベースの出現確率を推定したところ、 おおよそ4%程度であろうということがわかってきました。 購入者を分けずに全体のデータで推定した結果と大きくは変わらない結果となりました。

ということで、同僚Bが不正をしているという疑惑を持たずに済んでよかった。

【トップに戻る】

参考資料

  1. データ解析のための統計モデリング入門

  2. Pythonによるベイズ統計モデリング

    Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド

    Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド

    pymc3を使った入門書。 理論面も多少書かれており、こちらで手を動かしながら学習するのが良いと思う。

*1:最初の二つは、pymc2を使っているので情報が少し古いかも

*2:単純に楽しいので

*3:独立に確率が決まるという考えよりもこのような考え方の方が自然ですよね