2018-04-30

第５５回チョコボール計測

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

計測結果

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

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

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

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

基礎集計

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

項目	数
計測データ数	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%信用区間（上下それぞれ５％）は、下側が2.63%、上側が11.21%という結果です。
f:id:hippy-hikky:20180429102419p:plain

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

第５４回チョコボール計測

計測チョコボール

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

計測結果

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%信用区間（上下それぞれ５％）は、下側が1.12%、上側が8.06%という結果です。
f:id:hippy-hikky:20180429094822p:plain

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

第５３回チョコボール計測＋おまけ

計測チョコボール

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

計測結果

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%信用区間（上下それぞれ５％）は、下側が1.08%、上側が8.54%という結果です。
f:id:hippy-hikky:20180424234030p:plain

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

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

第５２回チョコボール計測＋おまけ

計測チョコボール

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

計測結果

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

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

基礎集計

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

項目	数
計測データ数	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%信用区間（上下それぞれ５％）は、下側が1.27%、上側が9.14%という結果です。
f:id:hippy-hikky:20180422235108p:plain

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

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

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

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

データ

同僚A、B、筆者の２人分の計測データ数です。カッコの中は銀のエンゼルの数です。
あのデータを示して以降、同僚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時点の予測結果。購入数が多い３人の結果だけに限定。

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

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

第５１回チョコボール計測

計測チョコボール

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

計測結果

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%信用区間（上下それぞれ５％）は、下側が1.27%、上側が9.86%という結果です。
f:id:hippy-hikky:20180422095533p:plain

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

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

チョコボール分析

概要
はじめに
データ
アプローチ
モデル
- チョコボールの個数の変化点予測
- 1個あたり平均重量の変化点予測
変化点の予測結果
- チョコボールの個数の変化点予測結果
- 1個あたり平均重量の変化点予測結果
まとめ
付録
参考文献
広告

概要

ここ最近、チョコボールの粒が大きくなったような気がする
ベイズ推定により、変化点の検出を試行
どうやら、賞味期限が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日

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

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

【トップに戻る】

データ

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

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

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

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

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

データは本来であれば、製造年月日でデータを並べるのが正しいのですが、製造日がわからないので、賞味期限でデータを並べています。（ただ、賞味期限は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個あたり平均重量の変化点予測

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

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

${ \displaystyle ｗ_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入門

作者: キャメロンデビッドソン=ピロン,玉木徹
出版社/メーカー: 森北出版
発売日: 2017/04/06
メディア: 単行本（ソフトカバー）
この商品を含むブログ (1件) を見る

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

【トップに戻る】

第５０回チョコボール計測

計測チョコボール

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

計測結果

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%信用区間（上下それぞれ５％）は、下側が1.43%、上側が10.45%という結果です。
f:id:hippy-hikky:20180417004122p:plain

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

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

計測結果

基礎集計

エンゼル出現確率の予測

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告

計測結果

基礎集計

エンゼル出現確率の予測

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

データ

出現確率予測結果

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告

概要

はじめに

データ

アプローチ

モデル

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

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

変化点の予測結果

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

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

まとめ

付録

参考文献

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告