【状態空間モデル】系列データとして銀のエンゼルの出現確率を推定する・改
【概要】
- 前回は、チョコボールの銀のエンゼル出現確率を時系列データのモデリングとして推論してみました
- しかしながら前回のモデルは、状態の依存関係が無いモデルになっていたので、今回はこの部分を改良してみました
- モデル改良の結果、ばらつきが少なく、直感にも合うような結果になりました
- ここ数ヶ月の間で確率が低下してる傾向が見えているような気のせいなような(要検証)
【目次】
はじめに
チョコボールのエンゼル出現確率を推論するにあたって、出現確率は一定(i.i.d)であると仮定するケースがほとんどだと思います*1。当ブログのメインコンテンツでもi.i.d仮定で出現確率を推論してきました(参考)。しかしながら、エンゼル出現確率は一定である保証はありません*2。
これに対して前回の記事では、エンゼル出現確率を系列データのモデリングとして拡張してみました。
しかし前回のモデルでは、状態の依存関係がなく、事前分布だけに系列としての情報を持っていました。これは、ベイズ的な逐次更新をパーティクルフィルタで実現していたとみなせるのかなと思います。
ということで今回は、「状態」の依存関係を明示的に持たせたモデルに改良してみます。
↓前回記事 chocolate-ball.hatenablog.com
モデリング
前回のモデルの課題
前回定式化したモデルを整理すると、以下の通りとなります。
(1)はシステムモデルで、時刻tの確率()は時刻t-1の時点の事後分布から生成されるとしました。(2)は観測モデルで、エンゼルの出現数は確率の二項分布に従うとしています。
この図のように、状態の時間的な依存性がなく、事前分布という形でのみ系列としての情報を伝播していくモデルになっていました。時間的に独立な試行に対して、逐次的に一定の確率を予測するという形(i.i.d仮定)になってしまっているので、「系列データ」としての予測モデルとは言えなかったのかなと思っています。
今回はこの部分を改良して、状態の時間的な依存性を考慮したモデリングを行います。
モデルの改良
今回扱うモデルは、上図のように状態空間モデルとしておなじみのモデルを考えます。
は時刻tでの状態(状態とは何かについては後述)、は時刻tでのエンゼル出現数です。前回同様に、チョコボールの仕様から時刻tは「賞味期限月」を表します。
前回との違いは「状態」が意味するものです。前回は、状態は出現確率を直接扱っていました。今回は状態と出現確率の関係は次のようになっています*3。
システムモデル()としては、前回同様に、出現確率は基本的には変化しないと考えます。そのため、状態はにガウスノイズで揺らぎを持たせることにします。
観測モデルは前回と同じで二項分布を利用します。
以上をまとめると、今回扱う状態空間モデルは次の通りです。
実装
推論には前回同様に逐次モンテカルロ法(パーティクルフィルタ)を利用します。利用したソースコードは私のgitで公開しています。リンクなどは都度貼るので必要な方は参考にしてもらえたらと思います(もし役に立ってくれたらうれしいです)。
観測モデル
観測モデルは二項分布です。なお二項分布のパラメータである確率は状態をシグモイド関数で変換しています。
def sigmoid(z): s = 1.0 / (1.0 + np.exp(-1. * z)) return s def obs_model(y_t, x_t, params=None): n_t = y_t[0] z = y_t[1] theta_t = sigmoid(x_t) w = stats.binom(n=n_t, p=theta_t).pmf(z) return w
システムモデル
状態はガウスノイズで揺らぎを持たせるとしますので、システムモデルは前回と比較してだいぶシンプルになりました。
def system_model(x_pre, params): scale = params['scale'] xt = stats.norm(loc=x_pre, scale=scale).rvs() return xt
パーティクルフィルタ
推論の計算は(たぶん)解析的な計算もそんなに苦労せずにできると思いますが、式の展開を考える手間を計算機に任せることにしました*4。
パーティクルフィルタの実装については、(ほぼ)前回の記事で書いた通りです。ちょっとだけ変更していまして、前回使ったノンパラメトリックな経験分布を使うケースと、今回のように既知の確率分布だけしか使わないケースに対応するようにしています。以下のリンクでソースコードを公開しています。
推論に利用した全コードは本記事末に貼っています。
実験
実際に計測結果を利用して実験してみます。
データ
2020/4/2現在の計測数と銀のエンゼルの出現数は次の通りです。
エンゼル | 個数 |
---|---|
なし | 693 |
銀 | 32 |
計 | 725 |
賞味期限月毎の単純なエンゼル出現割合を図にすると次の通りです。
薄目で見ると、2020年後半から2021年あたりから確率が下がったように見えなくもないです(気のせいかも)。でも、初期の購入数が多いため、単純には比較できなさそうです。
初期分布
状態の初期分布は、ほとんど何も考えずに適当な分散のガウスノイズから作りました。一様分布にするとシグモイド関数を掛けた際に確率が小さい値と大きい値に二極化してしまうので。
で、今回作った初期分布はこれです。
確率に変換すると次の通りです。もう少し分散を大きくした方が良かったかも。
推論結果
ということで推論してみました。
パーティクルの分布を重ねるとこんな感じ。
この結果を見ると、出現確率は2021年付近から下がったように見えます。が、気のせいかもしれません。最終時点の推論結果として、中央値は約3.0%でした。
今回のこの結果はあくまでも私が計測したデータに基づいたものであり、この推論結果が正しくエンゼル出現確率を推論できているかは全くの不明です。そもそも検証のステップを踏んでいないですし。
前回のモデルとの比較
前回のモデルと今回のモデルでの推論結果を比較してみます。
今回のモデルを利用した場合、出現確率が下がっているように推論されました。一方、前回のモデルでの推論結果は、5%を超える値になっており、データと合わせてみてもちょっと高すぎる推論になっているような気がします。どちらが尤もらしいかは未検証です。
次に、パーティクルの5%~95%領域を重ねたのが以下の図です。
この結果から、今回のモデルの方が推論結果のばらつきが小さくなっている様子が見えてきます。(パラメータに依存すると思いますが。)
モデルの評価をやっていないので感覚的な評価になりますが、今回の系列モデルの方が、直感に合うような結果になっていると思います。
まとめ
ということで今回は、前回に引き続き、チョコボールの銀のエンゼルの出現確率を状態空間モデルを利用して推論しました。
直感的な評価で言うと、今回のモデルの方がばらつきが小さく、データに適合するような推論結果になっているんじゃないかと思います。
で、この結果を見ると、賞味期限月が2021年に入ってから確率が下がっているようにも見えなくありません。しかしながら、この傾向が確かにいえそうなのか、統計的に意味があるくらいの差なのかについては未検証です。 この辺りの調査については引き続き行っていきます。
参考文献
[1], 樋口, 予測にいかす統計モデリングの基本, 講談社, 2011
予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで (KS理工学専門書)
- 作者:樋口 知之
- 発売日: 2011/04/07
- メディア: 単行本
[2], S. Thrun et.al. , Probabilistic Robotics, The MIT Press
Probabilistic Robotics (Intelligent Robotics and Autonomous Agents series)
- 作者:Thrun, Sebastian,Burgard, Wolfram,Fox, Dieter
- 発売日: 2005/08/19
- メディア: ハードカバー
[3], C.M. Bishop[著], 元田ら[訳], パターン認識と機械学習(上), シュプリンガー・ジャパン, 2007
- 作者:C.M. ビショップ
- 発売日: 2012/04/05
- メディア: 単行本(ソフトカバー)