銀のエンゼルと金のエンゼルの出現確率をベイズ推定する（金と銀を合わせて推定）

【概要】

これまで本ブログでは、金のエンゼルと銀のエンゼルの出現確率の推定を日々行っています
しかしこれまでの推定は、簡単のために金と銀のエンゼルを独立に推定してきました
そこで本記事では、金と銀のエンゼルを合わせて推定するため、多項分布を利用したモデルを構築します

【目次】

はじめに
アプローチ
- 前提と仮定
- モデル
データ
実験
終わりに
参考文献（お世話になった書籍）
広告

はじめに

当ブログでは、日々のチョコボール開封結果に基づいて金のエンゼルと銀のエンゼルの出現確率を推定しています。これまでの推定では、簡単のために、金のエンゼルと銀のエンゼルの出現確率をそれぞれベルヌーイ試行と仮定して独立に推定してきました（ベイズ推定でエンゼルの出現確率を予測する - チョコボール統計）。

しかしこのモデルは、現実を正しく表現していません。

チョコボールのエンゼルは金のエンゼルと銀のエンゼルが同時に現れることは無いですよね。「金のエンゼル」、「銀のエンゼル」、または、「ハズレ」の3状態のいずれかが実現します。そのため、金のエンゼルの出現確率と銀のエンゼルの出現確率を同時に推定する必要があります *1 *2。

また、「金のエンゼル2倍キャンペーン」という期間があったのですが、これも簡単のために、この期間のデータを除外して推定してきました。単純にデータが減るだけなので、この期間のデータも積極的に利用していきたいと思います。

ということで、これまで単純化するためにおいていた仮定を見直し*3、より現実に即したモデルでエンゼルの出現確率を推定していこうと思います。

【トップに戻る】

アプローチ

具体的に3つの状態(金のエンゼル、銀のエンゼル、ハズレ)を合わせて推定するモデルを考えていきます。モデルを考える前に、前提、仮定を明示しておきます。

前提と仮定

金のエンゼル、銀のエンゼルはランダムに発生し、データの全期間を通してエンゼルの出現確率は一定であるとします（金のエンゼル2倍キャンペーンは除く）。また、フレーバーによって出現確率は変わらないと仮定します。

次に、金のエンゼル2倍キャンペーン中は、金のエンゼルの出現確率が通常の出現確率の2倍であるとします（そのように公表されているので信じることにします）。

モデル

次に、チョコボールのエンゼルの発生過程をモデル化します。

まず、1章で述べたとおり実現値は3値｛金、銀、ハズレ｝です。そのため、多項分布を利用します。

${ \displaystyle m \sim {\rm Mult}(m | \pi, M) }$

ここで、 $m$ は全試行回数Mのうち、｛金、銀、ハズレ｝のそれぞれの個数を表すベクトルです。 $\pi$ は、｛金、銀、ハズレ｝の出現確率を表す3次元のベクトルで、このパラメータが推定対象です。

今回は、最尤推定ではなく、ベイズ推論を行いたいので、事前分布を設定する必要があります。そこで事前分布には、多項分布の共役事前分布であるディリクレ分布を利用します。

${ \displaystyle \pi \sim {\rm Dir}(\pi | \alpha) }$

ここで、 $\alpha$ はディリクレ分布のハイパーパラメータです。

ここまでを合わせて、Kruschkeのダイアグラム*4にすると以下のようなモデルとなります。

f:id:hippy-hikky:20190709230449p:plain — 図1．単純な多項分布のモデル．多項分布のパラメータの事前分布にはディリクレ分布を利用する．

次に、金のエンゼル2倍キャンペーンのデータを利用する方法を考えます。実はこのアイディアはネタが丸かぶりしている以下の記事を参考にさせていただきました（参考にというかほぼそのままです…）。

bayesmax.sblo.jp

上記の記事では、このキャンペーン期間のデータには確率に重みが付くというモデルにされています。それぞれの事象の重みを $w_*$ とすると、多項分布のパラメータ $\pi$ は以下のベクトルとなります。

${ \displaystyle \pi = (\pi_{none} w_{none}, \pi_{silver} w_{silver}, \pi_{gold} w_{gold})^{T} }$

ここで、重み $w_*$ は以下の値とします。

期間	$w_{none}$	$w_{silver}$	$w_{gold}$
通常期間	1.0	1.0	1.0
2倍キャンペーン	1.0	0.0	2.0

データ

今回利用するデータは、2017年11月~2019年7月までに当ブログ内で開封した566箱が対象です。なおグアムで購入したチョコボールについては、金のエンゼルも銀のエンゼルも共に存在しないため、対象外としています*5。

データをまとめると以下の通りです。 2.1説で説明した仮定により、推定対象のパラメータ（エンゼルの出現確率）は金のエンゼル2倍キャンペーン中の商品か否かにのみ依存するため、以下のように2つの期間に分けたデータとしました。

キャンペーン	ハズレ	銀のエンゼル	金のエンゼル
通常	432	20	0
金2倍	113	0	1

実験

パラメータ推論

2.2節に示した多項分布モデルのパラメータを推論します。 2.2節で述べたとおり今回の実験では、事前分布には共役事前分布であるディリクレ分布を利用します。そのため、ベイズの定理に従って事後分布を計算すると以下の通りディリクレ分布になります*6。

${ \displaystyle p(\pi | D ) \propto p(D | \pi)p(\pi) \\ \ = {\rm Mult}(m | \pi, M){\rm Dir}(\pi | \alpha) \\ \ = {\rm C} \prod^{K}_{k=1}\pi_k^{\alpha_k+m_k-1} = {\rm Dir}(\pi | \alpha+m) }$

ここで、 $\pi$ は事象の発生確率のベクトル（ここでは3次元ベクトル）、 mはデータを表し、各事象の発生回数を並べたベクトルで、 Mはデータの総数を表します( $\sum{m}=M$ )。 $\alpha$ はディリクレ事前分布のハイパーパラメータで、今回は適当な値を設定します。 $C$ は定数項を表します。

ということなのですが、今回はあえてPyMC3*7を利用し、サンプルによる近似事後分布を求めます（MCMC）*8。単純に私がPyMCを使いたかったのと、事前分布に共役ではない事前分布を設定できる柔軟さがあるので、今回は近似事後分布を求めました*9。

具体的なコードは、以下を参照ください。

実験結果

3章で示したデータを利用して、金のエンゼルと銀のエンゼルの出現確率を推定した結果を示します。 2章で述べたとおり、金のエンゼル2倍キャンペーンを含めないモデルと含めるモデルをそれぞれ推定しました。

金のエンゼル2倍キャンペーンを除いた場合

まず、問題を単純にするために金のエンゼル2倍キャンペーンを除いた場合の結果です。

f:id:hippy-hikky:20190709003916p:plain — 図x1. 各状態のパラメータの推定事後分布．左から，ハズレ，銀のエンゼル，金のエンゼルのそれぞれの出現確率事後分布

図x1は、多項分布のパラメータの事後分布をプロットしたものです。 3つの図はそれぞれ左から、ハズレ，銀のエンゼル，金のエンゼルの出現確率を表しています。それぞれの図には、期待値（mean）と95%HPD*10が記載されています。図x1から、95%の確率でハズレを引き、銀のエンゼルは約4.6%（2.8%~6.6%の間）で出現するであろうと推定できていることがわかります。一方、金のエンゼルは0.2%（0%~0.7%の間）であると推定しています*11。ここで利用したデータでは金のエンゼルは一つも出現していません。そのため最尤推定では0%と推定することになってしまいますが、事後分布では0.7%以下になるだろうと柔軟な予測ができています。

全てのデータを含めた場合

次に、2章で述べたように、エンゼルの出現確率に重みを載せて、全てのデータを利用して推定した結果を示します。

f:id:hippy-hikky:20190709005602p:plain — 図x2. 各状態のパラメータの事後分布(全データ利用)．左から，ハズレ，銀のエンゼル，金のエンゼルのそれぞれの出現確率事後分布

図x2は全てのデータを利用して、多項分布のパラメータの事後分布を推定した結果です。図のそれぞれの位置関係は、先の図x1と同じです。

図x1と図x2を比較すると、ほぼ同じような結果なのですが、金のエンゼルの事後分布を見ていただくと、分布の形状が微妙に異なっていることがわかると思います。今回のデータでは、金のエンゼル2倍キャンペーンの期間で金のエンゼルが一つ出ています*12。つまり確率0%では無いことははっきりしているため、0の付近の山が少し下がり、分布が右にシフトしている様子が確認できます。なお、図x1と同様に95%HPDの下限は0.0と表示されていますが、正確には、 $9.1e^{-3}\%$ と推定されました。

銀のエンゼルについては、4.2.1で利用したデータに加えてデータが追加されてはいますが、追加データについては銀のエンゼルの出現確率が0%として重みを設定しているため、推定結果には影響を与えていないということがわかります。

いくら買ったらエンゼルが当たるのか？

エンゼルの出現確率が推定できたのなら、次は、何個買えばエンゼルが当たるのかの見積もりが気になりますよね。

確率 $\theta$ のベルヌーイ試行において、n回成功するまでにk回失敗する確率を表現した分布として、「負の二項分布」が知られています (参考1, 参考2)。

${ \displaystyle f(k) = NB(k|\theta, n) = \binom{k+n-1}{n-1}\theta^{n}(1-\theta)^{k} }$

この分布を利用することで、エンゼルを得るまでに必要なチョコボールの個数を見積もることができそうです。しかし、上記の負の二項分布の定義を見ると、確率はpとして一点を代入する式になっていますが、 4.2節の推定で得られた結果は確率分布（事後分布）でした。事後分布で得られた推定結果の期待値*13を使って予測することもできますが、この確率分布にはデータがまだ十分でないための曖昧性が表現されているため、代表点で推定することは避けたいです。

そのため、以下の積分を計算することで、事後分布を利用した予測結果を得ることができます。

${ \displaystyle p(k | D) = \int p(k|\theta)p(\theta|D)d\theta \\ \ = \int NB(k|\theta, n)p(\theta | D) }$

$p(\theta | D)$ は4.2節で推定した事後分布です。期待値を計算するということですね。

ここで、今手元にある事後分布 $p(\theta | D)$ はサンプル集合として得られていることを思い出します。サンプル集合のためこのままでは上記の期待値を計算することはできません *14。しかし、サンプル集合で事後分布を予測できているため、サンプルごとの平均で積分を計算することができます。

${ \displaystyle p(k | D) \simeq \frac{1}{M}\sum^{M} NB(k|\theta_m, n) }$

ここで、Mはサンプルの数で、 $\theta_m$ はm番目の $\theta$ のサンプルを表します。

では早速金のエンゼル1枚と銀のエンゼルを5枚出すために必要なチョコボールの購入数を見積もって見ます。

銀のエンゼルを5つ得るまでに必要なチョコボールの購入数

f:id:hippy-hikky:20190709015215p:plain — 図x3．銀のエンゼルを5つ手に入れるまでに必要なチョコボールの個数の分布．

図x3は銀のエンゼルを5つ得るまでに必要な個数の分布（累積確率）です。事後分布を使って推定した結果（青線）と事後分布の期待値を使って推定した結果（赤線）を載せています。この図から、100個程度のチョコボールを買うことで、銀のエンゼルが5個得られる確率が50%を超えそうだということがわかります。また、銀のエンゼルの予測は、期待値を使った場合も事後分布を使った場合も概ね同じ程度であることがわかります。

金のエンゼルを5つ得るまでに必要なチョコボールの購入数

f:id:hippy-hikky:20190709015306p:plain — 図x4．金のエンゼルを1つ手に入れるまでに必要なチョコボールの個数の分布．

次に、図x4は金のエンゼルを1枚得るまでに必要な個数の分布（累積確率）です。こちらの図でも事後分布を使って推定した結果（青線）と事後分布の期待値を使って推定した結果（赤線）を載せています。この図から、金のエンゼルを得るためには、250個ほど買うことで50%を超えるということがわかります。 1,000個も買えば80%の確率で金のエンゼルが当たるという予想になっています。

期待値を使って予測した結果と事後分布を使って予測した結果を比較してみると、期待値を使って予測した方がポジティブな予測になっているのがわかります。図x2の事後分布を確認すると、金のエンゼルは右に裾が長い分布になっているため、期待値が少し高めなのだろうということがわかります。

終わりに

以上本記事は、金のエンゼルと銀のエンゼルを合わせて推定してみました。

結果としては、これまでの計測記事で示している独立に推定した場合とほぼ変わらないのですが、金のエンゼルは0.3%、銀のエンゼルは4.6%という推定値となりました。

また、エンゼルを得るまでに必要なチョコボールの購入数の見積もりを計算することで、事後分布を活用したメリットに触れることができました。

しかし、金のエンゼルはまだ一つで、しかも「金のエンゼル2倍キャンペーン」という特別な期間に出ただけなので、正確に推定するにはもっと多くのデータが必要です。そのため、今後も計測を続けていこうと思います。

参考文献（お世話になった書籍）

パターン認識と機械学習上
パターン認識と機械学習上
- 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
- 出版社/メーカー: 丸善出版
- 発売日: 2012/04/05
- メディア: 単行本（ソフトカバー）
- 購入: 6人クリック: 33回
- この商品を含むブログ (20件) を見る
ベイズ的な手法を軸にした機械学習の参考書。
Pythonによるベイズ統計モデリング
Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド
- 作者: オズワルドマーティン,Osvaldo Martin,金子武久
- 出版社/メーカー: 共立出版
- 発売日: 2018/06/22
- メディア: 単行本
- この商品を含むブログを見る
pymc3を使った入門書。理論面も多少書かれており、こちらで手を動かしながら学習するのが良いと思う。