系列データとして銀のエンゼルの出現確率を推定する（状態空間モデル）

【概要】

チョコボールの銀のエンゼル出現確率を時系列データとして扱い、状態空間モデルを使って推論してみた
状態空間モデルの推論はパーティクルフィルタを利用し、潜在変数の確率モデルとしてノンパラメトリックに推定した密度関数を利用してみた
結果としては毎月だいたい5%くらいの出現確率で一定していると予想できた

【目次】

はじめに
モデリング
- 前提/仮定
- モデル設計
  - ベータ分布の利用
  - 数値的な推論
実装
実験
- データ
- 推論結果
まとめ
参考文献

はじめに

当ブログではこれまで、チョコボールのエンゼル出現確率は、ある確率 $\theta$ で一定であると仮定してきました*1。そしてその仮定の元で、データから出現確率の予測を行っています。

chocolate-ball.hatenablog.com

しかし実際には、エンゼル出現確率が一定である保証はありません。実際、エンゼル「2倍キャンペーン」というキャンペーンが行われていたこともありました。

そこで、出現確率は時期によって変動する可能性を考慮し、時系列の状態推定問題として扱うことを考えてみます。 $\theta$ としていたのを $\theta_t$ として、時間による変化を許容するモデルに拡張するということです。

【トップに戻る】

モデリング

前提/仮定

【対象について】

推定の対象として「銀のエンゼル」の出現確率のみを扱います。またそれに伴って、「金のエンゼル2倍キャンペーン」期間のデータを除外します*2。

フレーバーについては、全ての味を同列で扱います。味毎に出現確率は変わらないということです。ですが、「パイナップル味」は除外します*3。これは国内向けの商品ではないので、エンゼルキャンペーンの対象外と明記されています。

【時間 $t$ について】

出現確率を時系列のモデリングで推定するとしましたが、エンゼルは「同一時期に製造された商品は同じ確率」と考えることにします。

「製造された時期」についてですが、チョコボールのパッケージには、製造日が記載されていません。賞味期限月が記されているのみであり、これが明確に時期を確認できる基準になります。購入日で分けても製造日と比例しているわけではありません。

モデル設計

状態空間モデル(参考文献[1, 2])は下図の通りです。

f:id:hippy-hikky:20210219222751p:plain

賞味期限月 $t$ の銀のエンゼルの出現確率を $\theta_t$ とし、エンゼルの出現数を $y_t$ とします。

銀のエンゼルの出現数 $y_t$ は確率 $\theta_t$ の二項分布に従うと仮定します。

ここで $n_t$ は賞味期限月 $t$ でのチョコボール購入数です。

次に、エンゼル出現確率 $\theta_t$ は基本的には変化しない（一定である）とします。

ここで、システムモデル（出現確率の遷移確率 $p(\theta_t | \theta_{t-1})$ ）の具体的な構造を設計します。二つのアイディアがあります。

観測モデルが二項分布なので、ベータ分布を使う
数値的に分布を定義する（ノンパラメトリックな推定密度関数を使う）

ベータ分布の利用

二項分布と共役関係にある分布としてベータ分布が知られています。

ここでも、出現確率のモデルとしてベータ分布を利用することで、解析的に推論結果を得ることができます（たぶん*4）。

エンゼル出現確率の推論モデルはとてもシンプルなモデルですので、通常は解析解を導出した方が良いと思います。

数値的な推論

ベータ分布のように既知の分布を仮定せず、KDEなどのノンパラメトリックな密度推定法(参考文献[3])を利用して数値的に確率分布を推定し、その推定した密度関数を利用します。既知の確率分布を仮定しないので、柔軟な分布構造を設計することができます。

数値的に密度関数が推定できれば、逆関数法を使ってサンプルを取得できます(参考文献[1])。

数値的な推論の場合、解析的に推論結果が得られずサンプルの集合として事後分布（フィルタ分布）を推論することになるので、パーティクルフィルタやアンサンブルカルマンフィルタなどを使う必要があります。*5

今回は、他への応用を考えて数値的に推論するこちらのアプローチを採用します*6。

なお利用する技術については、引用している参考文献を読んでいただくのが確実ですが、以下の記事でこれらの技術についてまとめていますので興味あれば参照いただけたらと思います。

learning-with-machine.hatenablog.com

【トップに戻る】

実装

状態空間モデルの状態推定には、逐次モンテカルロ法の数値的な実装であるパーティクルフィルタ（参考文献[2]）を利用します。パーティクルフィルタの実装については、別ブログの以下の記事で書いたクラスを（だいたいそのまま）利用します。

learning-with-machine.hatenablog.com

実装コードの全体はnotebookを添付しますのでよければ参照してもらえたらと思います。

観測モデルとシステムモデル

観測値（エンゼルの出現個数）は二項分布に従うと仮定しているので、二項分布の確率質量関数を利用するだけです。実装にはscipyを使いました。

システムモデルは前の時刻の事後分布をそのままシステムモデルにすると設計しました。今回は上述の通りパーティクルフィルタを利用して状態空間モデルの状態推定を行います。そのため、パーティクルの集合という形で事後分布が得られます。そこで、2.2節で書いた通り、数値的な手法を使って確率密度関数の推論とそこからのサンプルを行います。

パーティクルフィルタ

パーティクルフィルタについては、上述の通り、参考文献[1, 2]や私の書いたこちらのブログを参考にしていただけたらと思います。

ただし、以前私が書いた実装では、システムモデルと観測モデルは既知の確率分布に従うものとしていました。そのため、事後分布の数値的な密度推定を各時刻で行う必要があります。この部分を追加して次のように実装しました（長くなるのでリンクを貼っときました）。

particle_filter.py · GitHub

実装コード全体

推論の全体は次のnotebookを参照ください。

【トップに戻る】

実験

データ

今回利用するデータは2017年11月から当ブログで不定期で計測しているデータです。前処理として以下の操作を行っています(2.1節参照)。

賞味期限月の記録が漏れているデータを除外
パイナップル味（エンゼルキャンペーンの対象外）を除外
金のエンゼル 2倍キャンペーンの商品を除外（銀のエンゼルが出ない）

データ数は次の通りです（計695個）。

エンゼル	個数
なし	664
銀のエンゼル	31

賞味期限月毎の個数と各賞味期限月のエンゼルの出現割合は以下の通りです。

f:id:hippy-hikky:20210220225332p:plain — 賞味期限月毎の個数とエンゼル出現割合．オレンジの棒は各月の購入数，青線が各月のエンゼル出現割合．途切れているところはデータが欠損している月．

推論結果

上記のデータを利用して出現確率 $\theta_t$ を推論します。

推論結果は次の通りでした。

f:id:hippy-hikky:20210220225626p:plain — エンゼル出現確率の推論結果として、パーティクル集合の中央値(赤線)を追記．

この結果から、だいたい5%から7%程度であると推論していることがわかります。

次に、パーティクルの集合を重ねてみます。

f:id:hippy-hikky:20210220225840p:plain — 緑の点でパーティクル集合を追記．

わかりにくいですね。。。遠目でみると、データが欠損しているところではパーティクルが広がっているなーということは見えそうです。また、だいたい5%付近にパーティクルが集まっている（確率分布のピークがある）といえそうです（言えるか？）。

また、10ヶ月目くらいまではパーティクルが大きく拡がっており、ここまでの推論結果は信頼出来なさそうです。いわゆるburn in期間ということですね。

Chocolate Ball Viewerや最近の計測記事をみると、全体データをi.i.d仮定で予測した場合には、平均が4.6%程度と推論しています。なので、今回の推論結果はちょっと高めになっているように見えます。これが何に起因しているのかについては考察ができていません。（何か気付いた方、優しく教えてください）

【トップに戻る】

まとめ

ということで今回は、チョコボールの銀のエンゼル出現確率を系列データとして扱い、状態空間モデルをパーティクルフィルタを利用して推論してみました。また、潜在変数の確率モデルとしてノンパラメトリックに推定した密度関数を利用してみました。

推論結果としては、だいたい毎月5%程度で一定しているのかなと思います。一部10%くらいまでブレている月もありますが、ここはノイズ的なものと見れるのかなと考えています。

今回は一方通行での推論しか行っていませんので、平滑化を行うことでもう少し安定した推論結果になるのかなと思いますが、それは今後の課題とします。また、i.i.d.として推論した結果よりも少し高めの推論結果になっているように見えますので、この辺りも今後モデルの設計や実装にミスがないか確認していきたいと思います。

【トップに戻る】