チョコボール統計

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

第263回 チョコボール計測(ピーナツ)

こんばんは、チョコボール統計研究所です。
今夜もピーナツ味を3箱計測します。

f:id:hippy-hikky:20210224163056j:plain

こちらも更新してます → Chocolate Ball Viewer

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2021-02-24 2021-05-01 33.317 4.635 16 小山工場 BOX なし 28.682 1.793
2021-02-24 2021-05-01 33.08 4.641 17 小山工場 BOX なし 28.439 1.673
2021-02-24 2021-05-01 33.634 4.625 16 小山工場 BOX なし 29.009 1.813

今日もエンゼルさんは現れませんでした。

f:id:hippy-hikky:20210224163207j:plain

基礎集計

ピーナツ味の集計です。

項目
計測データ数 589
銀のエンゼル出現数 22
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.032 29.310 32.232 29.266
個数 14.000 16.000 20.000 16.312

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

この図は正味の重量のヒストグラムです。 赤い縦線が仕様(28g)を表しています。 青い太線で正規分布と仮定した最尤推定量をプロットしています。

エンゼル出現確率の予測

これまでに得られているデータを利用して金と銀のエンゼルの確率を推定します。 推定方法の詳細は以下の記事を参照ください。

銀のエンゼルと金のエンゼルの出現確率をベイズ推定する(金と銀を合わせて推定) - チョコボール統計

これまでに取得したデータは次の通りです。

項目
計測データ数 812
銀のエンゼル出現数 31
金のエンゼル出現数 1

この結果を使ってベイズ推定によるエンゼルの出現確率推定を行います。

ハズレ、銀のエンゼル、金のエンゼルの確率の推定結果は以下のとおりです。

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

銀のエンゼルは、3.2%〜6.1%の間と推定しており、期待値は4.6%です。

金のエンゼルは、0.00%~0.49%の間と推定しており、期待値は0.21%です。

広告

Amazonの欲しいものリスト作ってみました。 チョコボールのカンパ募集中です。
チョコボールをカンパする

森永製菓 チョコボール<ピーナッツ> 28g×20箱

森永製菓 チョコボール<ピーナッツ> 28g×20箱

  • 発売日: 2016/03/01
  • メディア: 食品&飲料

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

【概要】

  • チョコボール銀のエンゼル出現確率を時系列データとして扱い、状態空間モデルを使って推論してみた
  • 状態空間モデルの推論はパーティクルフィルタを利用し、潜在変数の確率モデルとしてノンパラメトリックに推定した密度関数を利用してみた
  • 結果としては毎月だいたい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の二項分布に従うと仮定します。


\begin{align}
 y_t \sim \mathrm{Binom}(y | \theta_t, n_t)
\end{align}

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

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


\begin{align}
 p(\theta_t | \theta_{t-1}) = p(\theta_{t-1} | \theta_{t-2}, y_{1:t-1})
\end{align}

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

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

ベータ分布の利用

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

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

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

数値的な推論

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

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

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

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

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

learning-with-machine.hatenablog.com

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節参照)。

  1. 賞味期限月の記録が漏れているデータを除外
  2. パイナップル味(エンゼルキャンペーンの対象外)を除外
  3. 金のエンゼル 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.として推論した結果よりも少し高めの推論結果になっているように見えますので、この辺りも今後モデルの設計や実装にミスがないか確認していきたいと思います。

【トップに戻る】

参考文献

【トップに戻る】

*1:「エンゼル2倍キャンペーン」の期間を除いてですが

*2:エンゼル2倍キャンペーン期間は銀のエンゼル出ませんので

*3:グアムで買いました

*4:ちゃんと計算してないけどできるはず

*5:アンサンブルカルマンフィルタは正直よくわかってないです

*6:使いたいだけってやつです

第262回 チョコボール計測(ピーナツ)

こんばんは、チョコボール統計研究所です。
今夜もピーナツ味を3箱計測します。

f:id:hippy-hikky:20210218003837j:plain

こちらも更新してます → Chocolate Ball Viewer

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2021-02-17 2021-05-01 33.199 4.632 17 小山工場 BOX なし 28.567 1.680
2021-02-17 2021-05-01 32.958 4.638 17 小山工場 BOX なし 28.320 1.666
2021-02-17 2021-05-01 33.5 4.638 15 小山工場 BOX なし 28.862 1.924

今日はエンゼルさんは現れませんでした。

f:id:hippy-hikky:20210218003919j:plain

基礎集計

ピーナツ味の集計です。

項目
計測データ数 586
銀のエンゼル出現数 22
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.032 29.312 32.232 29.268
個数 14.000 16.000 20.000 16.312

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

この図は正味の重量のヒストグラムです。 赤い縦線が仕様(28g)を表しています。 青い太線で正規分布と仮定した最尤推定量をプロットしています。

エンゼル出現確率の予測

これまでに得られているデータを利用して金と銀のエンゼルの確率を推定します。 推定方法の詳細は以下の記事を参照ください。

銀のエンゼルと金のエンゼルの出現確率をベイズ推定する(金と銀を合わせて推定) - チョコボール統計

これまでに取得したデータは次の通りです。

項目
計測データ数 809
銀のエンゼル出現数 31
金のエンゼル出現数 1

この結果を使ってベイズ推定によるエンゼルの出現確率推定を行います。

ハズレ、銀のエンゼル、金のエンゼルの確率の推定結果は以下のとおりです。

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

銀のエンゼルは、3.2%〜6.1%の間と推定しており、期待値は4.6%です。

金のエンゼルは、0.00%~0.49%の間と推定しており、期待値は0.21%です。

広告

Amazonの欲しいものリスト作ってみました。 チョコボールのカンパ募集中です。
チョコボールをカンパする

森永製菓 チョコボール<ピーナッツ> 28g×20箱

森永製菓 チョコボール<ピーナッツ> 28g×20箱

  • 発売日: 2016/03/01
  • メディア: 食品&飲料

第261回 チョコボール計測(ピーナツ)

こんばんは、チョコボール統計研究所です。
今夜もピーナツ味を3箱計測します。

f:id:hippy-hikky:20210216160343j:plain

こちらも更新してます → Chocolate Ball Viewer

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2021-02-16 2021-05-01 33.178 4.624 15 小山工場 BOX なし 28.554 1.904
2021-02-16 2021-05-01 33.085 4.66 15 小山工場 BOX なし 28.425 1.895
2021-02-16 2021-05-01 32.96 4.614 16 小山工場 BOX なし 28.346 1.772

今日はエンゼルさんは現れませんでした。

f:id:hippy-hikky:20210216160521j:plain

基礎集計

ピーナツ味の集計です。

項目
計測データ数 583
銀のエンゼル出現数 22
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.032 29.312 32.232 29.272
個数 14.000 16.000 20.000 16.312

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

この図は正味の重量のヒストグラムです。 赤い縦線が仕様(28g)を表しています。 青い太線で正規分布と仮定した最尤推定量をプロットしています。

エンゼル出現確率の予測

これまでに得られているデータを利用して金と銀のエンゼルの確率を推定します。 推定方法の詳細は以下の記事を参照ください。

銀のエンゼルと金のエンゼルの出現確率をベイズ推定する(金と銀を合わせて推定) - チョコボール統計

これまでに取得したデータは次の通りです。

項目
計測データ数 806
銀のエンゼル出現数 31
金のエンゼル出現数 1

この結果を使ってベイズ推定によるエンゼルの出現確率推定を行います。

ハズレ、銀のエンゼル、金のエンゼルの確率の推定結果は以下のとおりです。

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

銀のエンゼルは、3.2%〜6.1%の間と推定しており、期待値は4.6%です。

金のエンゼルは、0.00%~0.49%の間と推定しており、期待値は0.21%です。

広告

Amazonの欲しいものリスト作ってみました。 チョコボールのカンパ募集中です。
チョコボールをカンパする

森永製菓 チョコボール<ピーナッツ> 28g×20箱

森永製菓 チョコボール<ピーナッツ> 28g×20箱

  • 発売日: 2016/03/01
  • メディア: 食品&飲料

第260回 チョコボール計測(ピーナツ)

こんばんは、チョコボール統計研究所です。
今夜もピーナツ味を3箱計測します。

f:id:hippy-hikky:20210215223431j:plain

こちらも更新してます → Chocolate Ball Viewer

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2021-02-15 2021-05-01 33.665 4.636 17 小山工場 BOX なし 29.029 1.708
2021-02-15 2021-05-01 33.186 4.634 16 小山工場 BOX なし 28.552 1.784
2021-02-15 2021-05-01 33.264 4.663 16 小山工場 BOX なし 28.601 1.788

今日はエンゼルさんは現れませんでした。

f:id:hippy-hikky:20210215223705j:plain

基礎集計

ピーナツ味の集計です。

項目
計測データ数 580
銀のエンゼル出現数 22
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.032 29.314 32.232 29.276
個数 14.000 16.000 20.000 16.317

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

この図は正味の重量のヒストグラムです。 赤い縦線が仕様(28g)を表しています。 青い太線で正規分布と仮定した最尤推定量をプロットしています。

エンゼル出現確率の予測

これまでに得られているデータを利用して金と銀のエンゼルの確率を推定します。 推定方法の詳細は以下の記事を参照ください。

銀のエンゼルと金のエンゼルの出現確率をベイズ推定する(金と銀を合わせて推定) - チョコボール統計

これまでに取得したデータは次の通りです。

項目
計測データ数 803
銀のエンゼル出現数 31
金のエンゼル出現数 1

この結果を使ってベイズ推定によるエンゼルの出現確率推定を行います。

ハズレ、銀のエンゼル、金のエンゼルの確率の推定結果は以下のとおりです。

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

銀のエンゼルは、3.2%〜6.1%の間と推定しており、期待値は4.6%です。

金のエンゼルは、0.00%~0.49%の間と推定しており、期待値は0.21%です。

広告

Amazonの欲しいものリスト作ってみました。 チョコボールのカンパ募集中です。
チョコボールをカンパする

森永製菓 チョコボール<ピーナッツ> 28g×20箱

森永製菓 チョコボール<ピーナッツ> 28g×20箱

  • 発売日: 2016/03/01
  • メディア: 食品&飲料

第259回 チョコボール計測(ピーナツ)

こんばんは、チョコボール統計研究所です。
今夜もピーナツ味を3箱計測します。

f:id:hippy-hikky:20210209220249j:plain

こちらも更新してます → Chocolate Ball Viewer

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2021-02-09 2021-05-01 33.113 4.615 16 小山工場 BOX 28.498 1.781
2021-02-09 2021-05-01 33.181 4.638 16 小山工場 BOX なし 28.543 1.784
2021-02-09 2021-05-01 33.233 4.643 16 小山工場 BOX なし 28.590 1.787

久しぶりの銀のエンゼルです!

f:id:hippy-hikky:20210209220452j:plain

基礎集計

ピーナツ味の集計です。

項目
計測データ数 577
銀のエンゼル出現数 22
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.032 29.316 32.232 29.279
個数 14.000 16.000 20.000 16.317

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

この図は正味の重量のヒストグラムです。 赤い縦線が仕様(28g)を表しています。 青い太線で正規分布と仮定した最尤推定量をプロットしています。

エンゼル出現確率の予測

これまでに得られているデータを利用して金と銀のエンゼルの確率を推定します。 推定方法の詳細は以下の記事を参照ください。

銀のエンゼルと金のエンゼルの出現確率をベイズ推定する(金と銀を合わせて推定) - チョコボール統計

これまでに取得したデータは次の通りです。

項目
計測データ数 800
銀のエンゼル出現数 31
金のエンゼル出現数 1

この結果を使ってベイズ推定によるエンゼルの出現確率推定を行います。

ハズレ、銀のエンゼル、金のエンゼルの確率の推定結果は以下のとおりです。

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

銀のエンゼルは、3.2%〜6.2%の間と推定しており、期待値は4.6%です。

金のエンゼルは、0.00%~0.50%の間と推定しており、期待値は0.22%です。

広告

Amazonの欲しいものリスト作ってみました。 チョコボールのカンパ募集中です。
チョコボールをカンパする

森永製菓 チョコボール<ピーナッツ> 28g×20箱

森永製菓 チョコボール<ピーナッツ> 28g×20箱

  • 発売日: 2016/03/01
  • メディア: 食品&飲料

第258回 チョコボール計測(ピーナツ)

こんばんは、チョコボール統計研究所です。
今夜はピーナツ味を3箱計測します。

f:id:hippy-hikky:20210208235617j:plain

こちらも更新してます → Chocolate Ball Viewer

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2021-02-08 2021-05-01 33.231 4.662 16 小山工場 BOX なし 28.569 1.786
2021-02-08 2021-05-01 33.518 4.667 16 小山工場 BOX なし 28.851 1.803
2021-02-08 2021-05-01 33.172 4.663 16 小山工場 BOX なし 28.509 1.782

エンゼルさん現れずです。

f:id:hippy-hikky:20210208235929j:plain

基礎集計

ピーナツ味の集計です。

項目
計測データ数 574
銀のエンゼル出現数 21
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.032 29.316 32.232 29.283
個数 14.000 16.000 20.000 16.319

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

この図は正味の重量のヒストグラムです。 赤い縦線が仕様(28g)を表しています。 青い太線で正規分布と仮定した最尤推定量をプロットしています。

エンゼル出現確率の予測

これまでに得られているデータを利用して金と銀のエンゼルの確率を推定します。 推定方法の詳細は以下の記事を参照ください。

銀のエンゼルと金のエンゼルの出現確率をベイズ推定する(金と銀を合わせて推定) - チョコボール統計

これまでに取得したデータは次の通りです。

項目
計測データ数 797
銀のエンゼル出現数 30
金のエンゼル出現数 1

この結果を使ってベイズ推定によるエンゼルの出現確率推定を行います。

ハズレ、銀のエンゼル、金のエンゼルの確率の推定結果は以下のとおりです。

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

銀のエンゼルは、3.0%〜6.0%の間と推定しており、期待値は4.5%です。

金のエンゼルは、0.00%~0.50%の間と推定しており、期待値は0.22%です。

広告

Amazonの欲しいものリスト作ってみました。 チョコボールのカンパ募集中です。
チョコボールをカンパする

森永製菓 チョコボール<ピーナッツ> 28g×20箱

森永製菓 チョコボール<ピーナッツ> 28g×20箱

  • 発売日: 2016/03/01
  • メディア: 食品&飲料