エンゼルを高確率で当てられる運の良さを評価(階層ベイズモデルによる個人差の推定)
【概要】
- 恒常的に「運」が良いなんて本当にあるのでしょうか?
- 銀のエンゼルの出現確率が、チョコボールの購入者によって差異が出るのかを調査することで、運の良さを定量的に評価してみます。
- 購入者毎に独立に出現確率を推定、一般化線形混合モデルを利用して”個人差”として「運」を推定、階層ベイズモデルで購入者毎の出現確率を推定と三つのモデルで評価してみます。
- 結果は、今回のデータでは「運」の大きさに統計的な差異があるとは言えないということがわかりました。
【目次】
はじめに
「あの人は運が良いよね」とか「俺は運が悪い」っていうセリフ、聞いたり言ったりすることよくありますよね?
でも、恒常的に「運」が良いなんて本当にあるのでしょうか?
本当はランダムな現象であるにも関わらず、認知バイアスにより成功(失敗)の記憶が強化されているだけではないでしょうか??
ということで今回は、 チョコボールの銀のエンゼルの出やすさ(出現確率)が購入者によって差があるのかを分析し、 購入者の運の強さを測れるのかを試してみます。
もし統計的に明らかに発生確率に差が生じているのなら、 それは本当に神秘的な何かか、不正が行われているのかのどちらかのはず。
(内容に間違いや気になった点がありましたら、ご指摘くださると助かります)
問題設定
当ブログでは、チョコボールを買ってきて計測することをライフワークとしています。 これまでにもちょくちょく言及していますが、 チョコボールの購入者は筆者のみではありません。 友人、同僚から寄付(チョコボール)をいただくことがよくあります。
その中でも、よく寄付してくれる同僚が二人います(同僚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%となっています。
購入者毎のエンゼル出現確率の推定
チョコボールはランダムに購入すると仮定し、銀のエンゼルが出るか出ないかはある一定の確率に従って 毎回ランダムに決まるとします(つまりベルヌーイ試行)。 そのため、基本的なモデルとして銀のエンゼルが出る確率をとした二項分布を考えます。
エンゼル出現確率の推定に利用したコードはひとまとめにして5.5に貼り付けておきます。 参考にしていただけたらと思います。
最尤推定量
まずは最尤推定量を確認しましょう。 最尤推定については、以下の記事などを参照ください。
チョコボールに学ぶ実践的ベイズ統計モデリング入門 #1 - チョコボール統計
上記の記事などを参照すると、二項分布の最尤推定量は頻度であることがわかるので、 4.1節の表に記載の通りです(図で示すと以下の通り)。
これだけ見ると、Bの運が良いように見えると思います。 しかし、この差は本当に差があるのでしょうか?偶然とは考えられないでしょうか?
これは点推定なので、この差がどのくらい意味があるかはまだわかりません。 検定などすれば差の評価はできますが、 以降の節ではベイズ推定を使ってエンゼルが出る確率の 事後分布を利用して評価していきます。
独立に出現確率を推定
ベイズ推定で(iは購入者を指すindex)を推定しましょう。
モデルは、購入者毎の購入数とエンゼル出現確率をパラメータに持った二項分布です。
ベイズ推定なので事前分布を仮定する必要があります。 事前分布にはベータ分布(下記)を採用します。
ベータ分布は二項分布の共役事前分布ですし、 との二つのパラメータで様々な形状を表現できます。
このモデルを図で表すと以下のようになります(手書きでごめんなさい)。
の事前分布であるベータ分布のパラメータを、 つまり無情報事前分布として推定すると、結果は以下のようになりました。
この結果から、MAP(事後確率最大の点)は5.1に記載の最尤推定とほぼ同じように見えますが、 それぞれの分布は重なっており、購入者毎の明確な差は見えなそうだといえます。
事後分布の差の分布も見てみましょう。
この図は、購入者二人ずつの事後分布の差の分布です。 0付近に分布が集中していれば、差はないと言えます。 この図から、B対Cで少し差の分布が0から離れていますが、95%信用区間の中に収まっていますので、 差があるとは言えないということがわかります。
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ブロックで囲んだ中でモデルを定義します。 上から、の事前分布として無情報事前分布を設定し、 二行目で尤度として二項分布を用いることにしています。 三行目で、乱数のシードを固定してMCMCサンプル5000取得しています。
8~10セルまでは、MCMCサンプルを可視化して事後分布を確認しています。 この中では、10セル目で、購入者毎の差の分布を可視化していますが、 事後分布をMCMCサンプルで得ているので、単純に差をとった分布を確認すれば良いというのが お手軽で良いですね。
購入者の個人差をパラメータに加えて推定
5.2では、購入者独立にエンゼル出現確率を推定しました。 5.2の結果だけでも、統計的に明確な差は見えないということがわかりましたが、 もう少しモデルを複雑化してみましょう*2。
エンゼルの出現確率は、森永さんの設定により、ある決まった値になっているものと仮定します(これをベースの出現確率と呼ぶ)。 もし人によって「運」の量のような物があるのなら、 ベースの出現確率に個人毎が持っている運の要素が加えられて調整されていると考えます。
これをモデル化してみましょう。
は確率なので、[0, 1]の範囲に収まる必要があります。 そのため、ロジット関数を使ってエンゼル出現の対数オッズがで表されると仮定します。 ここで、は森永さんが設定するベースの出現確率を決めるパラメータで、 は購入者iの運の強さを表します。
とに適当な事前分布(裾が広い一様分布状の正規分布)を仮定して推定した結果は以下の通りです。
この結果から、運の要素というのは非常に小さいということがわかりますし、購入者毎の差もほとんどなさそうなことがわかります。 の事後分布の差の分布を以下に示します。
B対Cを見ると多少差が大きいようですが、信用区間の中に収まっており、明確に差があるとは言えないようです。
森永さんの設定する真の出現確率については、以下のように予測できました。
この図の横軸は、対数オッズです。 確率に変換するには、シグモイド関数を通す必要があります。 結果として、95%信用区間は以下の表の通りとわかりました。
2.5%点 | 平均 | 97.5%点 |
---|---|---|
0.25% | 4.19% | 34.00% |
まだ範囲が広いですね。データをもっと集める必要があります。
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)
ロジット関数を展開すると、二項分布のパラメータはシグモイド関数で表現できるので、
二項分布にはp=pm.math.sigmoid(p+u)
という形で入力しています。
推定対象のパラメータはの3種類ですので、
それぞれの事前分布を設定しています。
なお、suは個人差を決めるパラメータの事前分布のパラメータです。
セル番号16,17では、対数オッズで予測されるパラメータを確率に変換するためにシグモイド関数を定義して、 ベースの出現確率を算出しています。
階層モデルとして予測
もう少し見方を変えて、別のモデルを立ててみましょう。 5.3節では、購入者毎のエンゼルの出現確率は、 ベースの出現確率に個人差が加わったものとしてモデルを作りました。
ここでは、森永さんの設定するベースの出現確率の確率分布から購入者毎の出現確率がサンプルされると考えます。 つまり、5.2節のように購入者毎に独立に出現確率が決まるのではなく、 ベースの確率になんらかの調整が加わって購入者毎の確率が決まると考えます*3。
これは、階層ベイズモデルを考えることでモデル化できます。
モデルを以下の図に示します。
階層事前分布のそれぞれの事前分布は裾の広い正規分布(無情報事前分布)として、 上記のモデルに従って購入者毎のエンゼル出現確率を推定した結果を以下に示します。
購入者毎の差の分布は次の通りです。
これらの図から、やはり購入者によって差はほとんど無いということがわかりました。
では、の事前分布、つまり、ベースの出現確率の分布を確認してみます。
この結果から、ベースの出現確率の期待値は4%ということがわかりました。
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では、の事前分布はと固定されていましたが、 ここではも推定対象です。 そのため、これらのパラメータの事前分布を定義しています。
セル番号21では、の事前分布をプロットするコードを書いています。 これは、参考文献2からほぼそのままコピーしてきました。
ソースコード
まとめ
以上、購入者毎のエンゼル出現確率に差があるのかを3つのモデルを使って推定しました。
最尤推定(5.1)だけでは、意味のあるほどの差なのかわからなかったのですが、 ベイズ推定により事後分布を推定(5.1, 5.2, 5.3)したところ、 購入者毎にエンゼルの出現確率が明らかに異なるとは言えないということがわかりました。
また、購入者毎の効果(運)を明確に分離するモデル(5.3と5.4)を使って森永さんが設定しているであろう 銀のエンゼルのベースの出現確率を推定したところ、 おおよそ4%程度であろうということがわかってきました。 購入者を分けずに全体のデータで推定した結果と大きくは変わらない結果となりました。
ということで、同僚Bが不正をしているという疑惑を持たずに済んでよかった。
参考資料
データ解析のための統計モデリング入門
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人 クリック: 163回
- この商品を含むブログ (29件) を見る
- pymc3を使った入門書。 理論面も多少書かれており、こちらで手を動かしながら学習するのが良いと思う。
Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド
- 作者: オズワルドマーティン,Osvaldo Martin,金子武久
- 出版社/メーカー: 共立出版
- 発売日: 2018/06/22
- メディア: 単行本
- この商品を含むブログを見る