チョコボール統計

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

第55回 チョコボール計測

本日の計測結果です。
近所のコンビニで4箱買ってきました。

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2018-04-29 2018-11-01 34.179 4.841 18 小山工場 コンビニ(さいたま市 なし 29.338 1.630
2018-04-29 2018-11-01 34.876 4.863 18 小山工場 コンビニ(さいたま市 30.013 1.667
2018-04-29 2018-11-01 34.380 4.836 17 小山工場 コンビニ(さいたま市 29.544 1.738
2018-04-29 2018-11-01 34.106 4.886 17 小山工場 コンビニ(さいたま市 なし 29.220 1.719

久しぶりに銀のエンゼルが出ました!!

まさか一気に2つ出てくるとは。。。

ちなみに、新企画のために箱を開いて撮影してます。 新企画の詳細は別の記事にて公表予定です。

基礎集計

この集計はピーナツ味のチョコボールの集計結果です。

項目
計測データ数 153
銀のエンゼル出現数 2
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.913 29.384 30.657 29.450
個数 14.000 16.000 20.000 16.451

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

エンゼル出現確率の予測

通常版のエンゼルの予測を行っていきます。 これまでの通常版パッケージの開封結果は次の通りです。

項目
計測データ数 77
銀のエンゼル出現数 4
金のエンゼル出現数 0

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間(上下それぞれ5%)は、下側が2.63%、上側が11.21%という結果です。
f:id:hippy-hikky:20180429102419p:plain

次に金のエンゼルの出現確率の推定です。 90%信用区間(上下それぞれ5%)は、下側が0.07%、上側が3.83%という予測になっています。
f:id:hippy-hikky:20180429102503p:plain

広告

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

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

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

第54回 チョコボール計測

本日の計測結果です。
会社近くのコンビニで2箱買ってきました。

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2018-04-26 2018-12-01 34.094 4.794 16 小山工場 コンビニ(千代田区 なし 29.300 1.831
2018-04-26 2018-12-01 33.716 4.790 15 小山工場 コンビニ(千代田区 なし 28.926 1.928

今回もエンゼルは現れませんでした!

基礎集計

この集計はピーナツ味のチョコボールの集計結果です。

項目
計測データ数 149
銀のエンゼル出現数 0
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.913 29.384 30.657 29.447
個数 14.000 16.000 20.000 16.423

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

エンゼル出現確率の予測

通常版のエンゼルの予測を行っていきます。 これまでの通常版パッケージの開封結果は次の通りです。

項目
計測データ数 73
銀のエンゼル出現数 2
金のエンゼル出現数 0

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間(上下それぞれ5%)は、下側が1.12%、上側が8.06%という結果です。
f:id:hippy-hikky:20180429094822p:plain

次に金のエンゼルの出現確率の推定です。 90%信用区間(上下それぞれ5%)は、下側が0.07%、上側が3.75%という予測になっています。
f:id:hippy-hikky:20180429094851p:plain

広告

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

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

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

第53回 チョコボール計測+おまけ

本日の計測結果です。
同僚から5箱もらったものを開封します。

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2018-04-23 2018-11-01 34.131 4.877 14 小山工場 ドラックストア(横浜市 なし 29.254 2.090
2018-04-23 2018-11-01 34.067 4.887 15 小山工場 ドラックストア(横浜市 なし 29.180 1.945
2018-04-23 2018-11-01 34.349 4.871 14 小山工場 ドラックストア(横浜市 なし 29.478 2.106
2018-04-23 2018-10-01 34.324 4.780 15 小山工場 ドラックストア(横浜市 なし 29.544 1.970
2018-04-23 2018-10-01 34.000 4.733 15 小山工場 ドラックストア(横浜市 なし 29.267 1.951

今回もエンゼルは現れませんでした!

基礎集計

この集計はピーナツ味のチョコボールの集計結果です。

項目
計測データ数 147
銀のエンゼル出現数 0
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.913 29.385 30.657 29.452
個数 14.000 16.000 20.000 16.435

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

エンゼル出現確率の予測

通常版のエンゼルの予測を行っていきます。 これまでの通常版パッケージの開封結果は次の通りです。

項目
計測データ数 71
銀のエンゼル出現数 2
金のエンゼル出現数 0

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間(上下それぞれ5%)は、下側が1.08%、上側が8.54%という結果です。
f:id:hippy-hikky:20180424234030p:plain

次に金のエンゼルの出現確率の推定です。 90%信用区間(上下それぞれ5%)は、下側が0.08%、上側が4.01%という予測になっています。
f:id:hippy-hikky:20180424234105p:plain

銀のエンゼルの出現確率予測分布がだいぶ収束してきましたね。
現在のところ期待値として約4%くらいの出現確率であろうと予測しています。

広告

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

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

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

第52回 チョコボール計測+おまけ

本日の計測結果です。
金曜日に同僚から2箱もらったものと私が2つ買ったものを開封します。

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2018-04-22 2018-11-01 34.336 4.837 17 小山工場 コンビニ(さいたま市) なし 29.499 1.735
2018-04-22 2018-11-01 34.707 4.839 19 小山工場 コンビニ(さいたま市) なし 29.868 1.572
2018-04-22 2018-12-01 33.999 4.811 15 小山工場 コンビニ(千代田区 なし 29.188 1.946
2018-04-22 2018-12-01 34.142 4.854 15 小山工場 コンビニ(千代田区 なし 29.288 1.953

エンゼルは現れませんでした。
金のエンゼル2倍キャンペーンが終わって、銀復活したはずなんですが、 本当に復活しているのか疑ってしまうくらい出ないですね〜(48箱連続ハズレ)。

基礎集計

この集計はピーナツ味のチョコボールの集計結果です。

項目
計測データ数 142
銀のエンゼル出現数 0
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.913 29.386 30.657 29.456
個数 14.000 16.000 20.000 16.500

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

エンゼル出現確率の予測

通常版のエンゼルの予測を行っていきます。 これまでの通常版パッケージの開封結果は次の通りです。

項目
計測データ数 66
銀のエンゼル出現数 2
金のエンゼル出現数 0

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間(上下それぞれ5%)は、下側が1.27%、上側が9.14%という結果です。
f:id:hippy-hikky:20180422235108p:plain

次に金のエンゼルの出現確率の推定です。 90%信用区間(上下それぞれ5%)は、下側が0.06%、上側が4.67%という予測になっています。
f:id:hippy-hikky:20180422235139p:plain

おまけ:購入者毎の予測確率

先日の記事(↓)にて、購入者毎の運の強さを予測してみました。 chocolate-ball.hatenablog.com

あの計測以後、同僚AとBから追加してもらったので、改めて最近の状況を見てみます。

データ

同僚A、B、筆者の2人分の計測データ数です。カッコの中は銀のエンゼルの数です。
あのデータを示して以降、同僚Aにかなり協力もらっています。 が、未だにエンゼル現れず!

購入者 2018-04-10 2018-04-22
同僚A 8(0) 20(0)
同僚B 6(1) 9(1)
筆者 22(1) 32(1)

出現確率予測結果

手法やコードについては上記の記事を参照ください。

最初に、4/10時点の結果を貼っておきます。 同僚Aの運のなさが現れていましたね。

f:id:hippy-hikky:20180409000900p:plain
2018-04-10時点での予測結果

次に、最近の結果を載せます。

f:id:hippy-hikky:20180423000308p:plain
2018-04-22時点の予測結果。購入数が多い3人の結果だけに限定。

同僚Aの予測分布について、0の付近への偏りが凄い!
ということで、同僚Aの運の無さが際立ってきたという結果でした。

ここまでくると、同僚Aには是非とも当ててほしくないものです。

広告

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

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

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

第51回 チョコボール計測

本日の計測結果です。
同僚二人から3箱ずつもらいました。

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2018-04-18 2018-10-01 34.783 4.805 15 小山工場 ドラックストア(横浜市 なし 29.978 1.999
2018-04-18 2018-10-01 34.081 4.773 15 小山工場 ドラックストア(横浜市 なし 29.308 1.954
2018-04-18 2018-10-01 34.048 4.760 15 小山工場 ドラックストア(横浜市 なし 29.288 1.953
2018-04-18 2018-11-01 34.264 4.815 16 小山工場 コンビニ(千代田区 なし 29.449 1.841
2018-04-18 2018-11-01 34.280 4.838 16 小山工場 コンビニ(千代田区 なし 29.442 1.840
2018-04-18 2018-11-01 35.486 4.829 16 小山工場 コンビニ(千代田区 なし 30.657 1.916

エンゼルは現れませんでした。
全然出ないですねー。。。

基礎集計

この集計はピーナツ味のチョコボールの集計結果です。

項目
計測データ数 138
銀のエンゼル出現数 0
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.913 29.386 30.657 29.456
個数 14.000 16.000 20.000 16.500

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

エンゼル出現確率の予測

通常版のエンゼルの予測を行っていきます。 これまでの通常版パッケージの開封結果は次の通りです。

項目
計測データ数 62
銀のエンゼル出現数 2
金のエンゼル出現数 0

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間(上下それぞれ5%)は、下側が1.27%、上側が9.86%という結果です。
f:id:hippy-hikky:20180422095533p:plain

次に金のエンゼルの出現確率の推定です。 90%信用区間(上下それぞれ5%)は、下側が0.09%、上側が4.33%という予測になっています。
f:id:hippy-hikky:20180422095602p:plain

広告

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

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

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

チョコボールの大きさの変化点を捉える

概要

  • ここ最近、チョコボールの粒が大きくなったような気がする
  • ベイズ推定により、変化点の検出を試行
  • どうやら、賞味期限が2018年6月から2018年7月に変わるタイミング(計測時期的には2017年12月中旬)のあたりで変化したらしいと予測できた
  • 個数が1~2個減り、チョコボールの粒の平均重量が約0.2g(11%)重くなっているらしい

はじめに

去年の11月から半年弱、チョコボールの開封を不定期にやっているのですが、 最近のチョコボールは粒が揃っていて、大きい気がします。
例えば、このように。

f:id:hippy-hikky:20180415220437j:plain:w200
2017年11月26日計測
f:id:hippy-hikky:20180415220453j:plain:w300
2017年4月14日

普通はこのように写真を並べて目視で評価するかもしれませんが、 目視では確認できる量に限界があります。 また、主観に左右されるので信頼できないですよね。 (上の写真じゃ撮影位置が違うし全然変化がわからない。。。)
そこで、差があるのか、いつから差が出たのかを統計的に評価してみます。

この問題は、時系列データの変化点検出に応用できます。

【トップに戻る】

データ

今回使うデータは、ピーナツ味のデータのみを使います。 他の味は機械的に作るのでバラつきがほとんど無いでしょうし、 開封している数もピーナツ味が圧倒的に多く、他は少ないので。

ピーナツ味の全体の集計は下表の通りです。

項目
データ数 129箱開封
平均内包個数 16.5個
1個あたりの平均重量の平均 1.79g

データの散布図を以下に示します。

f:id:hippy-hikky:20180415232110p:plain
左:全体重量の散布図。中央:個数の散布図。右:1個あたり平均重量の散布図。横軸は全て開封数。

データは本来であれば、製造年月日でデータを並べるのが正しいのですが、 製造日がわからないので、賞味期限でデータを並べています。 (ただ、賞味期限は1月単位でしか表示がないので、その間の並びは不定です)

上図を見ると、全体重量(左図)は変化していないのに、 チョコボールの数(中央図)は60箱目くらいで少なくなっているように見えます。 この傾向は1個あたり平均重量(右図)が顕著で、 一箱あたりのチョコボール数が減っているので、1粒あたりの平均重量が重くなっている(=大きくなっている)ことが見えてきます。

では、実際に統計的に変化が起きていると言えるのか検証してみます。

【トップに戻る】

アプローチ

変化点の検出には、毎度おなじみのベイズ推定を用います。

変化が1回だと仮定し、変化前後での母分布のパラメータを推定します。 また、変化の起こった日(何個目のチョコボールから変化したのか)の予測も合わせて行います。 ちなみに、本当に変化が起きていない場合には、予測される変化点前後のパラメータには差異が出て来ないはずです。

実はこの問題、参考文献に載っている例とほぼほぼ同じです。 詳しくは書籍を読んでいただけたらと。

【トップに戻る】

モデル

変化点の有無はチョコボールの個数についての分析と、 1個あたり平均重量の分析の2種類の分析を行うことで評価してみます。

チョコボールの個数の変化点予測

はじめに1箱あたりに入っているチョコボールの数が変化しているのかを検証するためのモデルを立てます。

i箱目のチョコボールの個数N_iはパラメータ\lambdaポアソン分布に従うと仮定します。

{ \displaystyle
N_i \sim Poi(\lambda)
}

変化点\tauポアソン分布のパラメータが\lambda_1から\lambda_2に変化したとします。 ちなみに、ポアソン分布の期待値はパラメータ\lambdaになりますので、 \tau箱目から、チョコボール数の平均が\lambda_1から\lambda_2に変化したということをモデル化したものになります。

\lambdaの事前分布には、指数分布を用います。

{ \displaystyle
\lambda_i \sim Exp(\alpha)
}

ここで、\alphaは事前分布のパラメータであり、あまり過敏にならなくても良いものらしいです。 データ数が少なすぎなければ、事前分布がずれていてもリーズナブルな事後分布に収束してくれるはずです。

\tauも予測対象なので、事前分布を決める必要があります。 しかし、いつ変化したかという情報は何も持っていないので、一様分布を事前分布とします。

以上をまとめてpymcのコードにしたものは、付録を参照ください。

1個あたり平均重量の変化点予測

次に、一粒あたりの重量が変化しているのかを検証するためのモデルを立てます。 これは、1粒あたりの重量を予測する問題を応用します(過去の記事参照)。

i箱目の1個あたり平均重量w_iN(\mu, \sigma^{2})正規分布に従うと仮定します。

{ \displaystyle
w_i \sim N(\mu, \sigma^{2})
}

上記個数の変化点予測の場合と同様に、 変化点\tau正規分布のパラメータが(\mu_1, \sigma_1^{2})から(\mu_2, \sigma^{2}_{2})に変化したと仮定します。

今回は、 \muの事前分布には中心が観測データの平均、分散を100とした正規分布を使ってみます。 \sigma^{2}の事前分布には一様分布を使うことにします。 以前の記事と同じ設定です。

変化点\tauは上記と同じで一様分布を用います。

こちらのコードについても付録を参照ください。

【トップに戻る】

変化点の予測結果

変化点を予測した結果を見てみます。

チョコボールの個数の変化点予測結果

個数の変化点を予測した結果が以下の図です。

f:id:hippy-hikky:20180418000422p:plain
上:lamba_1の事後分布, 中央:lambda_2の事後分布, 下:tauの事後分布

\lambda_1\lambda_2の差異があるかを確認したところ、 95%以上の確信度で差異があると予測しています。 変化点の前後で、約17個から15~16個に減っているという予測結果になりました。 総重量が変化していないので、個数が減ったということは、粒が大きくなったということです。

変化点については、 60箱目くらいから80箱目くらいで変化していると予測しています。 先に示した散布図からも同様のことが読み取れますね。

しかし、変化点の分布は広く、 変化点を正確には予測できないようです。 実際、上記散布図の中央の個数の散布図を見てみると、 目で見てもあまりはっきりと変化点を読み取ることが難しいですね。

1個あたり平均重量の変化点予測結果

平均重量の変化点を予測した結果が以下の図です。 分散の事後分布については、ここでは掲示していません。

f:id:hippy-hikky:20180418000533p:plain
上:mu_1の事後分布, 中央:mu_2の事後分布, 下:tauの事後分布

結果は、個数の変化点予測と同じように、60箱目から80箱目付近でチョコボールの傾向に差が出ているということを示しています。 同じものを見方の異なるデータで表現したに過ぎないので、結果が同じになるのは当然ですね。
しかし、下段の\tauの予測分布を比較してみると、明らかにこの場合の\tauの分布が特定の値(60箱目、74箱目)に集中していることがわかります。 これは、60箱目か74箱目で変化しているということに強い確信を持っていることを示しています。 2山できているので、変化点が一つであるとの仮定が間違っていた可能性が高いです。

変化点前後の事後分布から、 1個あたり約0.2g重くなっているということがわかりました。 全体の1個あたり重量の平均が約1.8gなので、11%も重くなった(≒粒が大きくなった)ということが言えそうです。

【トップに戻る】

まとめ

ということで今回は、 チョコボールの大きさが変化したのか?ということについて分析してみました。

結果は、60箱目くらいを境にチョコボールの各粒の大きさが大きくなっているということがわかりました。 時期的には、賞味期限が2018年6月から2018年7月に変わるタイミングのようです。 ピーナツの生産地や時期的な違いによるのかもしれませんが、そこまではわかりません。

時期的な問題だと、1年の周期で変化の傾向が見えてくるかもしれないですね。

森永さん、公開しても良い範囲で教えてください!

【トップに戻る】

付録

今回分析に使ったコードを添付しておきます。 誰かの参考になれば幸いです。

【トップに戻る】

参考文献

Pythonで体験するベイズ推論 PyMCによるMCMC入門

Pythonで体験するベイズ推論 PyMCによるMCMC入門

pythonMCMC触ってみるなら。 例題が豊富で例題をトレースするだけでも面白いです。 ただし、本の趣旨の通り、pymcを触ってみることに興味のある人向け。 理論的な詳細は別の書籍を読むべき。

【トップに戻る】

広告

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

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

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

第50回 チョコボール計測

本日の計測結果です。
同僚から1箱もらいました。
新しくチョコボール集計に参加してくれた同僚です。

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2018-04-16 2018-10-01 34.088 4.778 15 小山工場 コンビニ(千代田区 なし 29.310 1.954

エンゼルは現れませんでした。
最初の一つでいきなりエンゼルが出るという奇跡は起きなかったですね。

基礎集計

この集計はピーナツ味のチョコボールの集計結果です。

項目
計測データ数 132
銀のエンゼル出現数 0
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.913 29.384 30.508 29.445
個数 14.000 16.000 20.000 16.545

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

エンゼル出現確率の予測

通常版のエンゼルの予測を行っていきます。 これまでの通常版パッケージの開封結果は次の通りです。

項目
計測データ数 56
銀のエンゼル出現数 2
金のエンゼル出現数 0

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間(上下それぞれ5%)は、下側が1.43%、上側が10.45%という結果です。
f:id:hippy-hikky:20180417004122p:plain

次に金のエンゼルの出現確率の推定です。 90%信用区間(上下それぞれ5%)は、下側が0.13%、上側が4.87%という予測になっています。
f:id:hippy-hikky:20180417004150p:plain

予測分布、だいぶ狭まって来ましたね。

広告

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

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

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