チョコボール統計

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

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

概要

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

第49回 チョコボール計測

本日の計測結果です。
同僚から2箱+自分で5箱買いました。

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2018-04-14 2018-11-01 34.376 4.817 16 小山工場 コンビニ(千代田区 なし 29.559 1.847
2018-04-14 2018-11-01 34.224 4.856 14 小山工場 コンビニ(千代田区 なし 29.368 2.098
2018-04-14 2018-09-01 33.955 4.723 15 小山工場 スーパー(さいたま市) なし 29.232 1.949
2018-04-14 2018-09-01 35.214 4.757 16 小山工場 スーパー(さいたま市) なし 30.457 1.904
2018-04-14 2018-09-01 34.000 4.728 16 小山工場 スーパー(さいたま市) なし 29.272 1.830
2018-04-14 2018-09-01 33.853 4.738 15 小山工場 スーパー(さいたま市) なし 29.115 1.941
2018-04-14 2018-09-01 34.060 4.738 15 小山工場 スーパー(さいたま市) なし 29.322 1.955

7箱開封するもエンゼルはれませんでした。

同僚からもらった分を撮影するの忘れてた。。。

基礎集計

この集計はピーナツ味のチョコボールの集計結果です。 今回、イチゴ味もあるんですが、イチゴ味の方は飛ばします。

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

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

エンゼル出現確率の予測

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

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

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

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

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

広告

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

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

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

第48回 チョコボール計測

本日の計測結果です。
同僚から2箱もらいました。

イチゴ味とピーナツ味の2箱です。

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2018-04-12 2018-11-01 34.046 4.824 15 小山工場 コンビニ(千代田区 なし 29.222 1.948
2018-04-12 2018-11-01 31.529 4.786 15 小山工場 コンビニ(千代田区 なし 26.743 1.783

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

今回協力してもらった同僚ですが、前回の分析記事の同僚Aです。 昨日も協力してもらったんですが、未だに連続ハズレ記録を更新中です。

基礎集計

この集計はピーナツ味のチョコボールの集計結果です。 今回、イチゴ味もあるんですが、イチゴ味の方は飛ばします。

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

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

エンゼル出現確率の予測

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

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

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

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

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

広告

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

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

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

第47回 チョコボール計測

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

パッケージがいつもと違うのに気づきましたか?
キョロちゃんセブンイレブンの制服着てます。

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

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2018-04-11 2018-10-01 34.526 4.789 16 小山工場 コンビニ(千代田区 なし 29.737 1.859
2018-04-11 2018-10-01 33.998 4.780 16 小山工場 コンビニ(千代田区 なし 29.218 1.826
2018-04-11 2018-10-01 33.889 4.748 16 小山工場 コンビニ(千代田区 なし 29.141 1.821

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

基礎集計

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

項目
計測データ数 123
銀のエンゼル出現数 0
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.913 29.389 30.508 29.446
個数 14.000 17.000 20.000 16.642

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

エンゼル出現確率の予測

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

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

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

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

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

広告

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

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

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

第46回 チョコボール計測

本日の計測結果です。
近所のコンビニで3箱買いました。

計測結果

date best_before weight box_weight number factory shop angel net_weight mean_weight
2018-04-10 2018-11-01 34.179 4.913 16 小山工場 コンビニ(さいたま市 なし 29.266 1.829
2018-04-10 2018-11-01 34.349 4.934 17 小山工場 コンビニ(さいたま市 なし 29.415 1.730
2018-04-10 2018-11-01 34.215 4.906 15 小山工場 コンビニ(さいたま市 なし 29.309 1.954

エンゼルは今日も現れませんでした。
金のエンゼル2倍キャンペーンが終わってからまだ一度も見ていないです。。。 そろそろ出て欲しい。。。

基礎集計

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

項目
計測データ数 120
銀のエンゼル出現数 0
金のエンゼル出現数 1
最小 中央値 最大値 平均
正味重量 28.913 29.389 30.508 29.448
個数 14.000 17.000 20.000 16.658

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

エンゼル出現確率の予測

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

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

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

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

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

広告

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

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

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

購入者によってエンゼルの出る確率は違う?

はじめに

運が強い弱いってこと、気になったことありますか?
チョコボールのエンゼルを含む「当たりクジ」をやたらと引いてしまう人がいるとかいないとか。 ゲン担ぎとか。

ということで今回は、チョコボールの購入者毎に、現時点での当たり確率を予測してみます。

【トップに戻る】

データ

データは、金のエンゼル2倍キャンペーンでは“無い”期間のチョコボールを対象にします。

ここまでの開封結果は下表の通りです。 同僚4人に協力してもらっています。

購入者 銀のエンゼル 総数
同僚A 0 8
同僚B 1 6
同僚C 0 2
同僚D 0 2
筆者 1 22

合計で40箱開封済みです。

頻度を見てみる

上記のデータを見てもらえばわかりますが、 筆者以外は購入数が10に満たない状態で、同僚A,C,Dはまだエンゼルが出ていません。

頻度を計算すると、同僚A,C,Dは0%、同僚Bは17%という値になります。

ここで、同僚AとC,Dでは購入数に違いがあります。 2箱しか開封していないC,Dに対して、Aは8箱も開封しています。
C,Dと比較してAの引きの弱さを感じないですか?

【トップに戻る】

少ないデータに対しての予測

本ブログでは、 このようにデータ数がまだ足りない場合の予測方法として、 「ベイズ推定」を使うアプローチをとってきました(以下の記事参照)。

chocolate-ball.hatenablog.com

ベイズ推定を使えば、 予測分布の広がり(分散)をもって、予測の信頼度を測ることができるのでした。
ということで今回も、ベイズ推定により購入者毎のエンゼル出現確率の予測をやってみます。

ベイズ推定による出現確率予測

ベイズ推定については、本ブログの他の分析記事や、参考文献等を参照ください。
上記のデータを使って、予測した結果が以下の図です。

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

この図は、以下の記事のコードをほぼ流用しています。 必要であればご参照ください。

ベイズ推定でエンゼルの出現確率を予測する - チョコボール統計

ベイズ推定で【銀の】エンゼルの出現確率を予測する - チョコボール統計

この図を見ると、 Aの予測分布が0に偏っているのに対し、 C,Dの分布がまだ広がっていることがわかります。 このことから、同僚Aの引きの弱さの直感が表現できていますね。

また、同僚Bの分布は広がってはいるものの、他と比較して高い位置に分布していることが見えます。

エンゼル出現数が0の同僚A,C,Dですが、 購入数の違いにより、Aの引きの弱さが分布に表現出来ているのが、ベイズ推定の面白いところですよね。

【トップに戻る】

まとめ

ということで、購入者別のエンゼル出現確率を予測してみました。

こういう分析は話のネタとしてはちょっとだけ笑えますが、 この結果はデータが少ないためのバラつきと考えられます。
結果の図を見てみると、 予測分布がいずれの購入者でもほぼ重なっている状態であり、 有意な差は見いだせません。

ランダムなサンプルの中からクジを引く場合、 当たりクジが入っている割合に依存して確率的に当たりが出るはずです。 (もしも確率に偏りが生じているなら、イカサマをしているという証拠)
しかし、今回のデータのようにサンプルが少ない場合には、 ばらつきのために運が良いとか悪いとか見えてしまう場合があります。 データを増やしていけば、この分布が収束していき、ほぼ重なるはずです。

もし重ならないようであれば、超常的な力が働いているのかも。。。

【トップに戻る】

おまけ

これだけだと、記事として軽いので、もう少しだけ分析を続けます。

データが少ない場合の予測は難しいものですが、 人間であれば常識や事前の知識に基づいてある程度のあたりをつけてデータを眺めることがあります。 良い場合と悪い場合あると思いますが、データが少ない場合には、大きく外す予測をしにくいという利点があります。

本ブログでは、 過去の記事で、 エンゼルの出現確率は不明なんだから無情報事前分布(一様分布)を使うということにしていました。

しかし今手元には、40箱開封して2個の銀のエンゼルが出たという証拠があります(5%)。

そこで、真の確率も5%から大きく外れた値ではないだろうというあたりをつけて予測をしてみたいと思います。

事前分布の設計

このような確率を扱う問題に対しての事前分布には、 ベータ分布が良く用いられます。 (ベータ分布については参考文献参照)

ベータ分布に従う確率変数は0から1の実数を取るので、確率をモデル化するのに調度良いという性質があります。

ベータ分布は、以下の式で表現され、パラメータは(\alpha, \beta)の2つです。

{ \displaystyle
f(x | \alpha, \beta) = \frac{x^{(\alpha - 1)}(1-x)^{(\beta-1)}}{B(\alpha, \beta)}
}

ここで、B(\alpha, \beta)はベータ関数を表します。

今回は、\alpha = 2.0\beta = 15.0として、約5%の値にピークを持つベータ分布を事前分布としてみました。(下図)

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

出現確率予測

上記のベータ分布を事前分布として、銀のエンゼルの出現確率の予測をしてみます。

実装は簡単で、過去の記事を基本にすると、 getMCMCResult関数の中で、事前分布の関数を以下のように変えているだけです。

# 出現確率pの事前分布
p = pm.Beta('p', alpha=2, beta=15)   ← この部分を変えているだけ
# 観測を結びつける
obs = pm.Bernoulli('obs', p, value=data, observed=True)

予測結果

予測結果は以下の図の通りです。

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

結論としては上記の結論とそう変わりはないのですが、 各購入者の予想分布が近づいています。
常識的に考えて、出現確率が50%を超えるようなことはまずありえないので、 妥当な結果なのかなと思います。

データが多量にあれば、事前分布の設計に気を付けることは必要ないのですが、 データが少なく、常識や事前知識がある程度ある場合には、 このように事前分布を適切なものに設計することで、 より早く妥当な分布に収束しやすくなると言えます。

また、ここでは言及しませんが、 今回のようにモデルにベルヌーイ分布を使っている場合、 ベータ分布は共役事前になり、事後分布もベータ分布になります。 つまりMCMCを使わずとも解析的に計算ができてしまいます。

このあたりの話も別の機会に書けたらと思います。

【トップに戻る】

参考文献

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

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

【トップに戻る】

広告

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

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

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