チョコボール統計

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

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

【概要】

  • 前回は、チョコボール銀のエンゼル出現確率を時系列データのモデリングとして推論してみました
  • しかしながら前回のモデルは、状態の依存関係が無いモデルになっていたので、今回はこの部分を改良してみました
  • モデル改良の結果、ばらつきが少なく、直感にも合うような結果になりました
    • ここ数ヶ月の間で確率が低下してる傾向が見えているような気のせいなような(要検証)

【目次】


はじめに

チョコボールのエンゼル出現確率を推論するにあたって、出現確率は一定(i.i.d)であると仮定するケースがほとんどだと思います*1。当ブログのメインコンテンツでもi.i.d仮定で出現確率を推論してきました(参考)。しかしながら、エンゼル出現確率は一定である保証はありません*2

これに対して前回の記事では、エンゼル出現確率を系列データのモデリングとして拡張してみました。

しかし前回のモデルでは、状態の依存関係がなく、事前分布だけに系列としての情報を持っていました。これは、ベイズ的な逐次更新をパーティクルフィルタで実現していたとみなせるのかなと思います。

ということで今回は、「状態」の依存関係を明示的に持たせたモデルに改良してみます。

↓前回記事 chocolate-ball.hatenablog.com

【トップに戻る】

モデリング

前回のモデルの課題

前回定式化したモデルを整理すると、以下の通りとなります。

\displaystyle{
\begin{align}
  \theta_t &\sim p(\theta_t | y_{1:t-1}) \cdots (1) \\
  y_t &\sim \mathrm{Binom}(y_t | \theta_t, n_t) \cdots(2)
\end{align}
}

(1)はシステムモデルで、時刻tの確率(\theta_t)は時刻t-1の時点の事後分布から生成されるとしました。(2)は観測モデルで、エンゼルの出現数y_tは確率\theta_tの二項分布に従うとしています。

f:id:hippy-hikky:20210405231955p:plain
前回のモデル.状態(θ)間に依存関係がない.

この図のように、状態\theta_tの時間的な依存性がなく、事前分布という形でのみ系列としての情報を伝播していくモデルになっていました。時間的に独立な試行に対して、逐次的に一定の確率\thetaを予測するという形(i.i.d仮定)になってしまっているので、「系列データ」としての予測モデルとは言えなかったのかなと思っています。

今回はこの部分を改良して、状態の時間的な依存性を考慮したモデリングを行います。

モデルの改良

f:id:hippy-hikky:20210405232753p:plain
状態空間モデル.今回は,状態の時間的な依存関係を持つように改良します.

今回扱うモデルは、上図のように状態空間モデルとしておなじみのモデルを考えます。

x_tは時刻tでの状態(状態とは何かについては後述)、y_tは時刻tでのエンゼル出現数です。前回同様に、チョコボールの仕様から時刻tは「賞味期限月」を表します。

前回との違いは「状態」が意味するものです。前回は、状態は出現確率を直接扱っていました。今回は状態x_tと出現確率\theta_tの関係は次のようになっています*3

\displaystyle{
\theta_t = \mathcal{f}(x_{t}) = \mathrm{Sigmoid}(x_t)
}

システムモデル(p(x_t | x_{t-1}))としては、前回同様に、出現確率は基本的には変化しないと考えます。そのため、状態x_tx_{t-1}ガウスノイズで揺らぎを持たせることにします。

\displaystyle{
x_t = x_{t-1} + \mathcal{N}(0, \sigma^2) \sim \mathcal{N}(x_t | x_{t-1}, \sigma^2)
}

観測モデルp(y_t | x_t)は前回と同じで二項分布を利用します。

以上をまとめると、今回扱う状態空間モデルは次の通りです。

\displaystyle{
\begin{align}
  x_t &\sim p(x_t | x_{t-1}) = \mathcal{N}(x_t | x_{t-1}, \sigma^2) \\\
  y_y &\sim p(y_t | x_t) = \mathrm{Binom}(y_t | \mathcal{f}(x_t), n_t)
\end{align}
}

【トップに戻る】

実装

推論には前回同様に逐次モンテカルロ法(パーティクルフィルタ)を利用します。利用したソースコードは私のgitで公開しています。リンクなどは都度貼るので必要な方は参考にしてもらえたらと思います(もし役に立ってくれたらうれしいです)。

観測モデル

観測モデルは二項分布です。なお二項分布のパラメータである確率は状態x_tシグモイド関数で変換しています。

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

パーティクルフィルタの実装については、(ほぼ)前回の記事で書いた通りです。ちょっとだけ変更していまして、前回使ったノンパラメトリックな経験分布を使うケースと、今回のように既知の確率分布だけしか使わないケースに対応するようにしています。以下のリンクでソースコードを公開しています。

particle_filter.py · GitHub

推論に利用した全コードは本記事末に貼っています。

【トップに戻る】

実験

実際に計測結果を利用して実験してみます。

データ

2020/4/2現在の計測数と銀のエンゼルの出現数は次の通りです。

エンゼル 個数
なし 693
32
725

賞味期限月毎の単純なエンゼル出現割合を図にすると次の通りです。

f:id:hippy-hikky:20210405235341p:plain
エンゼル出現割合(青線, 左軸)とチョコボール購入数(オレンジ棒, 右軸).2軸になっているので注意.

薄目で見ると、2020年後半から2021年あたりから確率が下がったように見えなくもないです(気のせいかも)。でも、初期の購入数が多いため、単純には比較できなさそうです。

初期分布

状態x_tの初期分布x_0は、ほとんど何も考えずに適当な分散のガウスノイズから作りました。一様分布にするとシグモイド関数を掛けた際に確率が小さい値と大きい値に二極化してしまうので。

で、今回作った初期分布はこれです。

f:id:hippy-hikky:20210405235953p:plain
状態xの初期分布.

確率\theta_0に変換すると次の通りです。もう少し分散を大きくした方が良かったかも。

f:id:hippy-hikky:20210406000017p:plain
状態xの初期分布にシグモイド関数をかけて確率に戻したもの.

推論結果

ということで推論してみました。

f:id:hippy-hikky:20210406000150p:plain
推論結果.先ほどのエンゼル出現割合に推論結果としてパーティクル集合の中央値を重ねた(赤線).

パーティクルの分布を重ねるとこんな感じ。

f:id:hippy-hikky:20210406000312p:plain
推論結果.観測値(エンゼルの出現割合, 青線),パーティクル集合(緑点),パーティクル中央値(赤線).

この結果を見ると、出現確率は2021年付近から下がったように見えます。が、気のせいかもしれません。最終時点の推論結果として、中央値は約3.0%でした。

今回のこの結果はあくまでも私が計測したデータに基づいたものであり、この推論結果が正しくエンゼル出現確率を推論できているかは全くの不明です。そもそも検証のステップを踏んでいないですし。

前回のモデルとの比較

前回のモデルと今回のモデルでの推論結果を比較してみます。

f:id:hippy-hikky:20210406001423p:plain
前回モデルとの比較.青線は各月のエンゼル出現割合.紫線が前回の独立なモデルでの推論結果,赤線が今回の系列モデルでの推論結果.推論結果はパーティクル集合の中央値.

今回のモデルを利用した場合、出現確率が下がっているように推論されました。一方、前回のモデルでの推論結果は、5%を超える値になっており、データと合わせてみてもちょっと高すぎる推論になっているような気がします。どちらが尤もらしいかは未検証です。

次に、パーティクルの5%~95%領域を重ねたのが以下の図です。

f:id:hippy-hikky:20210406001654p:plain
赤は今回の系列モデルでの推論結果.紫は前回のモデルでの推論結果.領域はパーティクル集合の5%から95%の範囲.

この結果から、今回のモデルの方が推論結果のばらつきが小さくなっている様子が見えてきます。(パラメータに依存すると思いますが。)

モデルの評価をやっていないので感覚的な評価になりますが、今回の系列モデルの方が、直感に合うような結果になっていると思います。

【トップに戻る】

まとめ

ということで今回は、前回に引き続き、チョコボール銀のエンゼルの出現確率を状態空間モデルを利用して推論しました。

直感的な評価で言うと、今回のモデルの方がばらつきが小さく、データに適合するような推論結果になっているんじゃないかと思います。

で、この結果を見ると、賞味期限月が2021年に入ってから確率が下がっているようにも見えなくありません。しかしながら、この傾向が確かにいえそうなのか、統計的に意味があるくらいの差なのかについては未検証です。 この辺りの調査については引き続き行っていきます。

【トップに戻る】

参考文献

【トップに戻る】

ソースコード全体

推論用notebook

パーティクルフィルタ

【トップに戻る】

*1:複数のyoutuberさんなどがスポットで調査した結果は私の知る限り全てこの仮定に則っています

*2:景品表示法?の関係で売り上げに対して景品の金額に上限が設けられているとか

*3:一般化線形モデルでいうところのロジスティック回帰、ロジットリンク関数みたいなものです

*4:こういうのが最近の流行りか?とも思ったけど単に僕が楽をしてるだけと言えます。