2019-03-25

第141回チョコボール計測（ピーナッツ）

本日の計測結果を報告します。今回も通常版のピーナッツを4箱計測します。
最近は1週間に一度のペースになってしまっているので、なんとか計測ペースを上げてデータを集めたいところです。

計測結果

date	best_before	weight	box_weight	number	factory	shop	angel	net_weight	mean_weight
2019-03-24	2019-08-01	34.001	4.720	15	小山工場	スーパー（さいたま市）	なし	29.281	1.952
2019-03-24	2019-08-01	34.816	4.748	17	小山工場	スーパー（さいたま市）	なし	30.068	1.769
2019-03-24	2019-08-01	33.862	4.754	16	小山工場	スーパー（さいたま市）	なし	29.108	1.819
2019-03-24	2019-08-01	34.001	4.723	16	小山工場	スーパー（さいたま市）	なし	29.278	1.830

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

基礎集計

通常版のピーナッツ味の集計結果です。

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

	最小	中央値	最大値	平均
正味重量	28.594	29.382	32.232	29.474
個数	14.000	17.000	20.000	16.549

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

エンゼル出現確率の予測

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

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

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間（上下それぞれ５％）は、下側が3.18%上側が6.68%という結果です。
f:id:hippy-hikky:20190325005514p:plain

現在の期待値は4.75%です。

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

今のところの金のエンゼルの期待値は0.25%です。

チョコボールに学ぶ実践的ベイズ統計モデリング入門 #2

チョコボール分析 Presentation

【概要】

3/16に表題のセミナーを開催しました
セミナーの資料の紹介と、資料に書ききれなかったことの補足をします

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

【目次】

はじめに
資料
セミナー概要
正味重量のモデル
最尤推定
ベイズ推定
階層ベイズモデル
まとめ
参考文献
広告

はじめに

2019/03/16に下記のセミナーを実施させていただきました。参加者のみなさん、ありがとうございました。

ml-for-experts.connpass.com

この記事では、セミナーのフォローアップを目的として、資料や資料の中で言い足りなかったことなどを補足します。

ちなみに、第一回は以下の通り開催しました。 ml-for-experts.connpass.com

第一回のフォローアップ記事は以下です。

チョコボールに学ぶ実践的ベイズ統計モデリング入門 #1 - チョコボール統計

【トップに戻る】

資料

当日使用した資料は下記Slideshareにアップしています。

20190316 bayes modeling from Yoichi Tokita

www.slideshare.net

ハンズオンで使用したコードは下記のGithub リポジトリに上がっております。 MITライセンスの範囲で自由に利用していただいて構いません。

github.com

【トップに戻る】

セミナー概要

実践的統計モデリングということで、前回はチョコボールの銀のエンゼルの出現確率についてモデリングを行いました*1。
第2回目の今回は、正味の重量の分布をモデリングしました。銀のエンゼルの出現有無は離散値（2値の確率変数）だったのですが、重量は連続値の確率変数になります。

前回と同様に、今回もまずは最尤推定で算出するとどのようになるのかを計算してみました。その後、ベイズ推定を行うことでどのように結果が得られるのかをハンズオンで体験してみました。

なおセミナー当日は、正味重量の推定に入る前に統計モデリングの目的について軽く復習しましたが、本記事ではその部分については触れません。気になる方は前回のフォローアップ記事を参照ください。

【トップに戻る】

正味重量のモデル

データ分析の基本は、まず可視化です。どのような問題*2かに依らず、取り組むデータがどのような構造をしているのかをはじめに体感しておくことは重要です。

と言うことで、今回のデータのヒストグラムは以下のようになりました。

f:id:hippy-hikky:20190321222812p:plain — チョコボールの正味重量のヒストグラム。左は正味の重量の分布でフレーバーによって基準の重さが異なるため分布の中央はそれぞれ異なる。右は、基準を統一するため、パッケージに表示の内容量からの差分をヒストグラムにしたもの。

今回は、簡単のために、重量分布は正規分布に従うと仮定しました。また実際の重量は、仕様上の重量よりも少し多めに入っていると仮定し、この量を「マージン」と呼ぶことにします。

${ \displaystyle p(y | \theta) = N(\mu, \sigma^{2}) \ }$

ここで、yは実際の重量を表す確率変数です。
平均 $\mu$ は仕様にマージンが加わったものとするので、

${ \displaystyle \mu = spec + \alpha }$

とします。 $spec$ はパッケージに表示の内容量で、仕様上の重量と呼びます。 $\alpha$ はマージンを表します。

しかし、可視化の結果を見ると正規分布とは言えなさそうです。セミナー内でも話しましたが、まずは最もシンプルなモデルから入り、不十分であればより複雑なモデルを構築していくことが良いと思います。

ということで以降では、正規分布をモデルとして考え、そのパラメータである平均と分散を推定していきます。パラメータの推定は第1回と同様に、「最尤推定」と「ベイズ推定」でそれぞれ推定してみます。

【トップに戻る】

最尤推定

最尤推定とは

最尤推定については、上記の資料や前回のフォローアップ記事、また、様々な書籍やブログ等でも紹介されています。

今回のモデルである正規分布のパラメータの最尤推定については、上記資料のp.25に記載しています。ざっくりと最尤推定について説明すると、尤度関数を最大化するパラメータを推定するという手法です。今回のモデルでは、データは独立同分布(i.i.d)であると仮定しているので、尤度関数はデータ毎の正規分布の同時確率になります（p.25参照）。尤度関数とは、データの当てはまりの良さを示す関数とイメージすると理解が進むかなと考えています。

最尤推定の計算

最尤推定の実際の計算も上記資料のp.25に記載しています。結果的には、標本平均と標本分散になります。

計算の細かいところについては、資料を参照ください。この計算は基礎なので、一度手で解かれることをお勧めします。

計算結果の詳細は、gitリポジトリにあるnotebookを参照ください。結果を図示すると以下のようになります。

f:id:hippy-hikky:20190322122059p:plain — 最尤推定量をデータに合わせて重ねた図。左図は重量の単純な推定分布（実線）。右図はマージンの分布。

ここまでの参考資料等

ここまでを理解するための資料としては、以下の書籍が参考になるかと思います。また、Google等で検索しても良い解説記事はかなりあるので、そちらを見ても良いと思います。

統計学入門

【トップに戻る】

ベイズ推定

モデルを正規分布と仮定したことで、最尤推定量は標本平均と標本分散であることが確認できました。最尤推定とは上記のように推定値を1点で表します（これを「点推定」と呼びます）。推定値なんだから、1点で解を得るのは当然じゃないかと思うかもしれませんが、点推定はいくつかの課題を抱えています。推定量が真のパラメータと一致するのはデータ数が無限大の極限であるので、データが少なければ偏りが生じます。

例えば、10個のサンプルデータの平均と1000個のサンプルデータの平均が一致したとしても、その値の信頼性は後者の方が高いですよね？点推定では、このように推定値の信頼性を認識できず、サンプルデータに過適合（オーバーフィッティング）してしまう恐れがあります。

ここに推定値の幅を合わせて推定できるベイズ推定を適用したいモチベーションがあると考えています。

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

ベイズ推定とは

ベイズ推定についての詳細は後述する参考資料を参照ください。

ざっくりと解説すると、以下のベイズの公式に基づいて、「事後分布」を推定するものです。

${ \displaystyle p(\theta | X) \propto p(X | \theta) p(\theta) }$

ここで、Xは計測データを表し、今回は $X=[x_1, x_2, \cdots ,x_n$ ]で各 $x_i$ は計測した重量を表します。右辺第１項は尤度関数であり、最尤推定で使ったものと同じものです。第２項がパラメータの事前分布で、これがパラメータの範囲を決めたり、事前の信念を与えるものです。

この式は、事後分布が尤度と事前分布の積に比例するという意味になります。最尤推定に対して、事前分布が正則化をしているという見方もできます。

ベイズ推定で正味重量の推定

実際に、正味重量のパラメータ推定をベイズ推定でやってみたいと思います。モデルの設定については資料のp.31,32に記載してあります。

モデルを考えていく際に、図を使って考えると理解が進み、改良点等も見えてきます。以下のような図を書きながら考えていきます(p.32)。

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

この図は、実際に計測する重量は正規分布からサンプルされるというモデルを示したもので、正規分布のパラメータ $\mu$ は $\mu=spec+\alpha$ であり、 $\alpha$ を実際には推定するものとしています。また、 $\alpha$ の事前分布には正規分布を利用することとしています。

このようにモデルを構築できると、あとはPPL(確率プログラミング言語)を利用すれば容易に解を得ることができます。今回のハンズオンでは、pymc3を利用しています。コードについては、githubのnotebookを参照ください。

ベイズ推定で正味重量の推定結果

ベイズ推定により正味重量を推定した結果、下図のようになりました。

f:id:hippy-hikky:20190322124755p:plain — ベイズ推定でマージンの量を推定した結果（事後分布）。上から、ピーナッツ、いちご、チョコバナナのマージン量の事後分布。

前節で記載のモデルの通り、実際の推定はマージン量を推定しています。この結果を見ると、ピーナッツとイチゴのマージン量に大きな差は無いようですが、チョコバナナは明らかにマージンが少ないと言えます。

今回は３つのフレーバーだけを推定しているのでわかりにくいのですが、他のフレーバーも推定してみると、最近発売の季節限定チョコボールのマージン量は今回のチョコバナナと同様の傾向があることがわかっています*3。しかし、ピーナッツ等レギュラー品には変化が無いようです。製造時のマージン量についてなんらかの方針転換があったのかなと推察していますが、実際のところはよくわかっていません。

ここまでの参考資料等

ハンズオンではPythonを利用し、Pythonの確率プログラミング言語であるpymc3を利用しました。 pymc3を利用した参考書としては、Pythonによるベイズ統計モデリングが参考になるかと思います。こちらは理論的な解説も載っているので、まずこちらの書籍を利用して触りながらベイズ統計を理解すると良いかもしれません。

ベイズ統計の実践的参考書としてはPythonで体験するベイズ推論もあり、こちらは事例が豊富で参考になるシーンは多いかと思います。しかし、こちらの書籍ではpymc2を利用しているため、コードの記述方法を多少読み替える必要があります。

【トップに戻る】

階層ベイズモデル

今回のセミナーでは最後に、上記のモデルを拡張する形で階層ベイズモデルについても少し触れました *4。

階層ベイズモデルとは

まず、6.2節に記載したベイズ推定のモデルの図を思い出します。これは、正規分布に従って重量がサンプルされるだろうとするモデルで、その正規分布の事前分布に正規分布を置くとしたモデルです。マージンを表すパラメータ $\alpha$ の事前分布には、平均 $\mu$ で分散が $\sigma^{2}$ の正規分布を仮定しています。

$\alpha$ についてまず注目します。 $\alpha_i$ はフレーバーiのマージン量を表しており、その事前分布に $N(\mu, \sigma^{2})$ を仮定しています。ここで、 $\alpha_i$ の事前分布 $N(\mu, \sigma^{2})$ には何の前提知識も持ち合わせていないので、ハンズオンでは分散の大きな正規分布を使って実行しました（Github上のnotebook参照）。

これは、 $\alpha_i$ には何の仮定も置かず、フレーバー毎のデータに基づいて独立に $\alpha_i$ を推定していることになります。事前分布は正則化の意味を持っているということを思い出すと、データが潤沢にある状況であれば問題は起こりませんが、データが少ないような場合には偏りが生じる可能性があります。

そこで、 $\alpha_i$ の事前分布もデータに基づいて値の範囲を狭めた強い仮定の事前分布にしたいという欲求が生まれます。このために、単純に前処理として、全データで平均と分散を推定するような方法も考えられますが、階層ベイズモデルにより一度に推定することができます。

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

上図は6.2節の図を拡張したもので、 $\alpha_i$ の事前分布にさらに事前分布が設定されていることがわかります。このような事前分布の事前分布を「ハイパー事前分布」や「階層事前分布」と呼びます。このハイパー事前分布の意味ですが、 $\alpha_i$ の共通した分布を表します。マージン $\alpha_i$ は何らかの形で製造装置等に設定するパラメータだろうと想像します。このパラメータがフレーバーと全く独立に決められているというのは考えにくく、会社としてのルールがあるのだろうと想像しますので、そのルールを推定するものであると見ることができます。

階層ベイズモデルというのは、このように事前分布にさらに事前分布を追加したモデルを指します。

階層ベイズモデルの実装

pymcのようなPPLを利用することで、階層モデルも容易に計算することができます。実際のコード例は、上記Githubの以下のnotebookを参照ください。

SeminarChocoBall/hierarchical_model.ipynb at master · tok41/SeminarChocoBall · GitHub

7.1節に図で書いたモデルが素直に記載されているのがわかるかと思います。

ここまでの参考資料等

階層ベイズモデルについては、先に紹介したPythonによるベイズ統計モデリングが参考になるかと思います。こちらはpymc3の実装について主に書かれているので、実際に手を動かしながら実感したい場合に有効です。

また、もう少し理論面を抑えたい場合には、データ解析のための統計モデリング入門もおすすめです。

なお、セミナーで紹介したnotebookのうちエンゼルの出現確率の推定については、当ブログの以下の記事に詳しく記載していますので合わせて参照いただけたらと思います。

chocolate-ball.hatenablog.com

【トップに戻る】

まとめ

「チョコボールに学ぶ実践的ベイズ統計モデリング入門」の第2回セミナーとして正味の重量の推定をやってみました。連続量の確率変数のモデリングの例でしたし、フレーバー毎のマージン量 $\alpha_i$ を推定するという複数のパラメータについてのモデリングの例でした。マージン量という比較できる量にして、その値の推定を行うことで、チョコバナナ味とその他の味の傾向が統計的に異なるということを推定しました。

このようなモデリングの応用例としては、装置の経時的な劣化を推定してアラートを出したり、異常検知のような問題に適用できる可能性があります。

【トップに戻る】

参考文献

統計学入門
統計学入門 (基礎統計学?)
- 作者: 東京大学教養学部統計学教室
- 出版社/メーカー: 東京大学出版会
- 発売日: 1991/07/09
- メディア: 単行本
- 購入: 158人クリック: 3,604回
- この商品を含むブログ (79件) を見る
有名な統計学の入門書。かっちり書いてあるので、統計学を学ぶにはすごく良い。
Pythonによるベイズ統計モデリング
Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド
- 作者: オズワルドマーティン,Osvaldo Martin,金子武久
- 出版社/メーカー: 共立出版
- 発売日: 2018/06/22
- メディア: 単行本
- この商品を含むブログを見る
pymc3を使った入門書。理論面も多少書かれており、こちらで手を動かしながら学習するのが良いと思う。
Pythonで体験するベイズ推論
Pythonで体験するベイズ推論 PyMCによるMCMC入門
- 作者: キャメロンデビッドソン=ピロン,玉木徹
- 出版社/メーカー: 森北出版
- 発売日: 2017/04/06
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (1件) を見る
pymcを利用した事例ベースの参考書。事例が豊富で参考になるシーンは出てくるかもしれないが、階層モデルについては記載がない。また、pymc2なので、pymc3と記法が少し違うので注意。
データ解析のための統計モデリング入門
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人クリック: 163回
- この商品を含むブログ (29件) を見る
一般化線形モデルを軸に階層ベイズモデルまで丁寧に記載されている。 Rとstanを利用した実装例もあるが、実装部分は無視しても理解できる。

【トップに戻る】

第140回チョコボール計測（ピーナッツ）

計測チョコボール

本日の計測結果を報告します。今回も通常版のピーナッツを4箱計測します。
一週間ほど開きましたが、ちょっといろいろ詰まっていたもので。

計測結果

date	best_before	weight	box_weight	number	factory	shop	angel	net_weight	mean_weight
2019-03-18	2019-10-01	34.002	4.706	17	小山工場	コンビニ（港区）	なし	29.296	1.723
2019-03-18	2019-10-01	34.080	4.719	16	小山工場	コンビニ（港区）	なし	29.361	1.835
2019-03-18	2019-09-01	34.181	4.740	15	小山工場	コンビニ（港区）	なし	29.441	1.963
2019-03-18	2019-09-01	33.682	4.671	16	小山工場	コンビニ（港区）	なし	29.011	1.813

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

基礎集計

通常版のピーナッツ味の集計結果です。

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

	最小	中央値	最大値	平均
正味重量	28.594	29.382	32.232	29.475
個数	14.000	17.000	20.000	16.555

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

エンゼル出現確率の予測

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

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

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間（上下それぞれ５％）は、下側が3.24%上側が7.94%という結果です。
f:id:hippy-hikky:20190318233436p:plain

現在の期待値は4.83%です。

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

今のところの金のエンゼルの期待値は0.27%です。

第139回チョコボール計測（ピーナッツ）

計測チョコボール

本日の計測結果を報告します。今回も通常版のピーナッツを4箱計測します。

引き続きこれです。

うちの妻から、バレンタインだってさ pic.twitter.com/E1ClZOuCee
— yoichi_t (@yoichi_t) 2019年2月14日

計測結果

date	best_before	weight	box_weight	number	factory	shop	angel	net_weight	mean_weight
2019-03-10	2019-08-01	33.998	4.735	15	小山工場	スーパー（さいたま市）	なし	29.263	1.951
2019-03-10	2019-08-01	34.084	4.756	14	小山工場	スーパー（さいたま市）	なし	29.328	2.095
2019-03-10	2019-08-01	34.242	4.708	15	小山工場	スーパー（さいたま市）	なし	29.534	1.969
2019-03-10	2019-08-01	34.304	4.757	15	小山工場	スーパー（さいたま市）	なし	29.547	1.970

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

基礎集計

通常版のピーナッツ味の集計結果です。

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

	最小	中央値	最大値	平均
正味重量	28.594	29.383	32.232	29.477
個数	14.000	17.000	20.000	16.561

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

エンゼル出現確率の予測

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

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

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間（上下それぞれ５％）は、下側が3.20%上側が6.84%という結果です。
f:id:hippy-hikky:20190310230046p:plain

現在の期待値は4.91%です。

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

今のところの金のエンゼルの期待値は0.25%です。

エンゼルを高確率で当てられる運の良さを評価（階層ベイズモデルによる個人差の推定）

チョコボール分析階層ベイズ GLMM

【概要】

恒常的に「運」が良いなんて本当にあるのでしょうか？
銀のエンゼルの出現確率が、チョコボールの購入者によって差異が出るのかを調査することで、運の良さを定量的に評価してみます。
購入者毎に独立に出現確率を推定、一般化線形混合モデルを利用して”個人差”として「運」を推定、階層ベイズモデルで購入者毎の出現確率を推定と三つのモデルで評価してみます。
結果は、今回のデータでは「運」の大きさに統計的な差異があるとは言えないということがわかりました。

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

【目次】

はじめに
問題設定
アプローチ
データ
- データの概要
購入者毎のエンゼル出現確率の推定
まとめ
参考資料

はじめに

「あの人は運が良いよね」とか「俺は運が悪い」っていうセリフ、聞いたり言ったりすることよくありますよね？
でも、恒常的に「運」が良いなんて本当にあるのでしょうか？本当はランダムな現象であるにも関わらず、認知バイアスにより成功（失敗）の記憶が強化されているだけではないでしょうか？？

ということで今回は、チョコボールの銀のエンゼルの出やすさ（出現確率）が購入者によって差があるのかを分析し、購入者の運の強さを測れるのかを試してみます。

もし統計的に明らかに発生確率に差が生じているのなら、それは本当に神秘的な何かか、不正が行われているのかのどちらかのはず。

（内容に間違いや気になった点がありましたら、ご指摘くださると助かります）

【トップに戻る】

問題設定

当ブログでは、チョコボールを買ってきて計測することをライフワークとしています。これまでにもちょくちょく言及していますが、チョコボールの購入者は筆者のみではありません。友人、同僚から寄付(チョコボール)をいただくことがよくあります。

その中でも、よく寄付してくれる同僚が二人います(同僚B,Cと呼ぶ、本記事の中でAは筆者を指す)。同僚Bは人格者で周りから好かれている良いやつです。一方、同僚Cはちょっとクセがあります。この二人、どちらも頻繁に寄付をしてくれるのですが、同僚Bからもらったチョコボールから銀のエンゼルがやたらと出る感じがします。（一方Cは。。。以下略）

そこで、筆者（A）、同僚B、同僚Cが購入したチョコボールのエンゼル出現確率に差があるのかを調査してみます。ランダムにチョコボールを購入しているのであれば、統計的に差異は生じないはずです。本当に差があるのだとしたら、それは神秘的な作用によるものか、または、不正（エンゼル出現に関するなんらかの情報を持っている）が行われているのか。。。

【トップに戻る】

アプローチ

チョコボールのエンゼル出現確率の推定は、ベイズ推定を利用します。エンゼル出現率とベイズ推定については当ブログの以下の記事を参照ください。*1

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

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

チョコボールに学ぶ実践的ベイズ統計モデリング入門 #1 - チョコボール統計

今回は、上記の記事のモデルをそのまま拡張して、購入者毎に独立に出現確率を推定してみます。その後モデルを拡張し、個人差を加えた階層モデルとしてエンゼルの出現確率を推定し、比較を行ってみます。それぞれのモデルの詳細は後述します。

【トップに戻る】

データ

利用するデータは当ブログで集計しているデータです。全データから、「金のエンゼル2倍キャンペーン」の期間のデータを除外します。このキャンペーンでは銀のエンゼルが含まれないと公式で発表されているからです。

データの概要

データの集計期間：2017/11/19〜2019/03/07

購入者	購入数	銀のエンゼル	頻度
A	142	6	0.042
B	80	5	0.063
C	59	2	0.034
total	281	13	0.046

全体での最尤推定量（頻度）は4.6%となっています。

【トップに戻る】

購入者毎のエンゼル出現確率の推定

チョコボールはランダムに購入すると仮定し、銀のエンゼルが出るか出ないかはある一定の確率に従って毎回ランダムに決まるとします（つまりベルヌーイ試行）。そのため、基本的なモデルとして銀のエンゼルが出る確率を $\theta$ とした二項分布を考えます。

${ \displaystyle p(y | \theta) = \binom{n}{y}\theta^{y}(1-\theta)^{n-y} }$

エンゼル出現確率の推定に利用したコードはひとまとめにして5.5に貼り付けておきます。参考にしていただけたらと思います。

最尤推定量

まずは最尤推定量を確認しましょう。最尤推定については、以下の記事などを参照ください。

チョコボールに学ぶ実践的ベイズ統計モデリング入門 #1 - チョコボール統計

データの性質を理解する - 機械と学習する

上記の記事などを参照すると、二項分布の最尤推定量は頻度であることがわかるので、 4.1節の表に記載の通りです（図で示すと以下の通り）。

f:id:hippy-hikky:20190309154627p:plain — 購入者A,B,C毎のエンゼル出現確率の最尤推定量

これだけ見ると、Bの運が良いように見えると思います。しかし、この差は本当に差があるのでしょうか？偶然とは考えられないでしょうか？

これは点推定なので、この差がどのくらい意味があるかはまだわかりません。検定などすれば差の評価はできますが、以降の節ではベイズ推定を使ってエンゼルが出る確率 $\theta$ の事後分布を利用して評価していきます。

ソースコード解説（クリックで展開）

最尤推定量の計算とプロットは、5.5のセル番号6です。
最尤推定量は頻度なので、単純にエンゼルの数を全体数で割っているだけなので、解説は不要かと思います。

独立に出現確率を推定

ベイズ推定で $\theta_i$ （iは購入者を指すindex）を推定しましょう。

モデルは、購入者毎の購入数 $n_i$ とエンゼル出現確率 $\theta_i$ をパラメータに持った二項分布です。

${ \displaystyle p(y_i | \theta_i) = \binom{n_i}{y_i}\theta_i^{y_i}(1-\theta_i)^{n_i-y_i} }$

ベイズ推定なので事前分布を仮定する必要があります。事前分布にはベータ分布(下記)を採用します。

${ \displaystyle p(\theta) = Beta(\alpha, \beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1} }$

ベータ分布は二項分布の共役事前分布ですし、 $\alpha$ と $\beta$ の二つのパラメータで様々な形状を表現できます。

このモデルを図で表すと以下のようになります（手書きでごめんなさい）。

f:id:hippy-hikky:20190309160816p:plain:w300 — 購入者毎に独立にエンゼル出現確率推定のモデル

$\theta$ の事前分布であるベータ分布のパラメータを $(\alpha, \beta)=(1,1)$ 、つまり無情報事前分布として推定すると、結果は以下のようになりました。

f:id:hippy-hikky:20190309161227p:plain — 購入者毎に独立に出現確率を推定した結果（ $\theta$ の事後分布）

この結果から、MAP（事後確率最大の点）は5.1に記載の最尤推定とほぼ同じように見えますが、それぞれの分布は重なっており、購入者毎の明確な差は見えなそうだといえます。

事後分布の差の分布も見てみましょう。

f:id:hippy-hikky:20190309161708p:plain — 事後分布の差の分布。左からA対B、A対C、B対C。

この図は、購入者二人ずつの事後分布の差の分布です。 0付近に分布が集中していれば、差はないと言えます。この図から、B対Cで少し差の分布が０から離れていますが、95%信用区間の中に収まっていますので、差があるとは言えないということがわかります。

ソースコード解説（クリックで展開）

購入者毎に独立にエンゼルの出現確率を推定している部分は、5.5のセル番号7です。

with pm.Model() as comparing_buyer_i:
    theta = pm.Beta('theta', alpha=1, beta=1, shape=len(set(buyer_idx)))
    angel = pm.Binomial('angel', n=total_counts[lst_buyer], p=theta[lst_buyer], observed=angel_counts[lst_buyer])
    trace_iso = pm.sample(5000, chains=1, random_seed=100)

pymc3の記法として、withブロックで囲んだ中でモデルを定義します。上から、 $\theta$ の事前分布として無情報事前分布を設定し、二行目で尤度として二項分布を用いることにしています。三行目で、乱数のシードを固定してMCMCサンプル5000取得しています。

8~10セルまでは、MCMCサンプルを可視化して事後分布を確認しています。この中では、10セル目で、購入者毎の差の分布を可視化していますが、事後分布をMCMCサンプルで得ているので、単純に差をとった分布を確認すれば良いというのがお手軽で良いですね。

購入者の個人差をパラメータに加えて推定

5.2では、購入者独立にエンゼル出現確率 $\theta_i$ を推定しました。 5.2の結果だけでも、統計的に明確な差は見えないということがわかりましたが、もう少しモデルを複雑化してみましょう*2。

エンゼルの出現確率は、森永さんの設定により、ある決まった値になっているものと仮定します（これをベースの出現確率と呼ぶ）。もし人によって「運」の量のような物があるのなら、ベースの出現確率に個人毎が持っている運の要素が加えられて調整されていると考えます。

これをモデル化してみましょう。

f:id:hippy-hikky:20190309163347p:plain:w300 — 購入者の個人差(運)を加えたモデル

$\theta$ は確率なので、[0, 1]の範囲に収まる必要があります。そのため、ロジット関数を使ってエンゼル出現の対数オッズが $p+u_i$ で表されると仮定します。ここで、 $p$ は森永さんが設定するベースの出現確率を決めるパラメータで、 $u_i$ は購入者iの運の強さを表します。

$p$ と $u_i$ に適当な事前分布（裾が広い一様分布状の正規分布）を仮定して推定した結果は以下の通りです。

f:id:hippy-hikky:20190309164118p:plain — 購入者毎の運を定めるパラメータ $u_i$ の事後分布

この結果から、運の要素というのは非常に小さいということがわかりますし、購入者毎の差もほとんどなさそうなことがわかります。 $u_i$ の事後分布の差の分布を以下に示します。

f:id:hippy-hikky:20190309164357p:plain — 運の要素の事後分布の差の分布。左からA対B、A対C、B対C。

B対Cを見ると多少差が大きいようですが、信用区間の中に収まっており、明確に差があるとは言えないようです。

森永さんの設定する真の出現確率については、以下のように予測できました。

f:id:hippy-hikky:20190309165721p:plain — 真のエンゼル出現確率(対数オッズ)の事後分布

この図の横軸は、対数オッズ $\log{\frac{q}{1-q}}$ です。確率に変換するには、シグモイド関数を通す必要があります。結果として、95%信用区間は以下の表の通りとわかりました。

2.5%点	平均	97.5%点
0.25%	4.19%	34.00%

まだ範囲が広いですね。データをもっと集める必要があります。

ソースコード解説（クリックで展開）

購入者毎の個人差を加えたモデルを計算している部分は、5.5のセル番号11です。

with pm.Model() as comparing_buyer_m1:
    su = pm.HalfNormal('su', sd=10)
    p = pm.Normal('p', mu=0, sd=10)
    u = pm.Normal('u', mu=0, sd=su, shape=len(set(buyer_idx)))
    
    angel = pm.Binomial('angel', 
                        n=total_counts[lst_buyer], 
                        p=pm.math.sigmoid(p+u[lst_buyer]), observed=angel_counts[lst_buyer])
    
    trace_h1 = pm.sample(3000, chains=3, random_seed=100)

ロジット関数を展開すると、二項分布のパラメータ $\theta$ はシグモイド関数で表現できるので、二項分布にはp=pm.math.sigmoid(p+u)という形で入力しています。推定対象のパラメータは $p, u_i, su$ の３種類ですので、それぞれの事前分布を設定しています。なお、suは個人差を決めるパラメータ $u_i$ の事前分布のパラメータです。

セル番号16,17では、対数オッズで予測されるパラメータを確率に変換するためにシグモイド関数を定義して、ベースの出現確率を算出しています。

階層モデルとして予測

もう少し見方を変えて、別のモデルを立ててみましょう。 5.3節では、購入者毎のエンゼルの出現確率は、ベースの出現確率 $p$ に個人差 $u_i$ が加わったものとしてモデルを作りました。

ここでは、森永さんの設定するベースの出現確率の確率分布から購入者毎の出現確率がサンプルされると考えます。つまり、5.2節のように購入者毎に独立に出現確率が決まるのではなく、ベースの確率になんらかの調整が加わって購入者毎の確率が決まると考えます*3。

これは、階層ベイズモデルを考えることでモデル化できます。

モデルを以下の図に示します。

f:id:hippy-hikky:20190309173044p:plain:w250 — 購入者毎のエンゼル出現確率は、真の確率(事前分布)からサンプルされると仮定するモデル

階層事前分布 $Beta(\alpha, \beta)$ のそれぞれの事前分布は裾の広い正規分布（無情報事前分布）として、上記のモデルに従って購入者毎のエンゼル出現確率を推定した結果を以下に示します。

f:id:hippy-hikky:20190309173720p:plain — 階層モデルで推定した購入者毎のエンゼル出現確率の事後分布

購入者毎の差の分布は次の通りです。

f:id:hippy-hikky:20190309173848p:plain — 階層モデルで推定した購入者毎のエンゼル出現確率の事後分布の差の分布

これらの図から、やはり購入者によって差はほとんど無いということがわかりました。

では、 $\theta$ の事前分布、つまり、ベースの出現確率の分布を確認してみます。

f:id:hippy-hikky:20190310220856p:plain — $\theta$ の事前分布。濃い青線が期待値を示す。オレンジはMCMCサンプルからのランダムなサンプル。

この結果から、ベースの出現確率の期待値は4%ということがわかりました。

ソースコード解説（クリックで展開）

購入者毎に独立ではなく、階層ベイズモデルでエンゼルの出現確率を推定している部分は、 5.5のセル番号18です。

with pm.Model() as comparing_buyer_m2:
    alpha = pm.HalfCauchy('alpha', beta=10)
    beta = pm.HalfCauchy('beta', beta=10)
    
    theta = pm.Beta('theta', alpha=alpha, beta=beta, shape=len(set(buyer_idx)))
    
    angel = pm.Binomial('angel', n=total_counts[lst_buyer], p=theta[lst_buyer], observed=angel_counts[lst_buyer])
    
    trace_h2 = pm.sample(2000, chains=3, random_seed=100)

5.2の購入者毎に独立に出現確率を推定したコードと似ています。 5.2では、 $\theta$ の事前分布は $Beta(\alpha=1,\beta=1)$ と固定されていましたが、ここでは $\alpha, \beta$ も推定対象です。そのため、これらのパラメータの事前分布を定義しています。

セル番号21では、 $\theta_i$ の事前分布をプロットするコードを書いています。これは、参考文献2からほぼそのままコピーしてきました。

ソースコード

まとめ

以上、購入者毎のエンゼル出現確率に差があるのかを3つのモデルを使って推定しました。

最尤推定(5.1)だけでは、意味のあるほどの差なのかわからなかったのですが、ベイズ推定により事後分布を推定(5.1, 5.2, 5.3)したところ、購入者毎にエンゼルの出現確率が明らかに異なるとは言えないということがわかりました。

また、購入者毎の効果（運）を明確に分離するモデル（5.3と5.4）を使って森永さんが設定しているであろう銀のエンゼルのベースの出現確率を推定したところ、おおよそ4%程度であろうということがわかってきました。購入者を分けずに全体のデータで推定した結果と大きくは変わらない結果となりました。

ということで、同僚Bが不正をしているという疑惑を持たずに済んでよかった。

【トップに戻る】

参考資料

データ解析のための統計モデリング入門
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人クリック: 163回
- この商品を含むブログ (29件) を見る
Pythonによるベイズ統計モデリング
Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド
- 作者: オズワルドマーティン,Osvaldo Martin,金子武久
- 出版社/メーカー: 共立出版
- 発売日: 2018/06/22
- メディア: 単行本
- この商品を含むブログを見る
pymc3を使った入門書。理論面も多少書かれており、こちらで手を動かしながら学習するのが良いと思う。

*1:最初の二つは、pymc2を使っているので情報が少し古いかも

*2:単純に楽しいので

*3:独立に確率が決まるという考えよりもこのような考え方の方が自然ですよね

2019-03-08

第138回チョコボール計測（ピーナッツ）

計測チョコボール

本日の計測結果を報告します。今回も通常版のピーナッツを4箱計測します。

引き続きこれです。

うちの妻から、バレンタインだってさ pic.twitter.com/E1ClZOuCee
— yoichi_t (@yoichi_t) 2019年2月14日

計測結果

date	best_before	weight	box_weight	number	factory	shop	angel	net_weight	mean_weight
2019-03-07	2019-08-01	34.114	4.777	14	小山工場	スーパー（さいたま市）	なし	29.337	2.096
2019-03-07	2019-08-01	34.600	4.762	15	小山工場	スーパー（さいたま市）	なし	29.838	1.989
2019-03-07	2019-08-01	34.044	4.743	15	小山工場	スーパー（さいたま市）	なし	29.301	1.953
2019-03-07	2019-08-01	33.952	4.756	16	小山工場	スーパー（さいたま市）	なし	29.195	1.825

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

基礎集計

通常版のピーナッツ味の集計結果です。

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

	最小	中央値	最大値	平均
正味重量	28.594	29.383	32.232	29.477
個数	14.000	17.000	20.000	16.581

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

エンゼル出現確率の予測

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

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

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間（上下それぞれ５％）は、下側が3.39%上側が7.03%という結果です。
f:id:hippy-hikky:20190308011240p:plain

現在の期待値は4.98%です。

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

今のところの金のエンゼルの期待値は0.26%です。

第137回チョコボール計測（ピーナッツ）

計測チョコボール

本日の計測結果を報告します。今日は通常版のピーナッツを4箱計測します。

これ↓です。

うちの妻から、バレンタインだってさ pic.twitter.com/E1ClZOuCee
— yoichi_t (@yoichi_t) 2019年2月14日

計測結果

date	best_before	weight	box_weight	number	factory	shop	angel	net_weight	mean_weight
2019-03-06	2019-08-01	34.107	4.747	15	小山工場	スーパー（さいたま市）	なし	29.360	1.957
2019-03-06	2019-08-01	33.999	4.753	16	小山工場	スーパー（さいたま市）	なし	29.246	1.828
2019-03-06	2019-08-01	34.134	4.752	15	小山工場	スーパー（さいたま市）	なし	29.382	1.959
2019-03-06	2019-08-01	34.001	4.736	14	小山工場	スーパー（さいたま市）	なし	29.265	2.090

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

基礎集計

通常版のピーナッツ味の集計結果です。

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

	最小	中央値	最大値	平均
正味重量	28.594	29.384	32.232	29.478
個数	14.000	17.000	20.000	16.599

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

エンゼル出現確率の予測

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

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

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

はじめに銀のエンゼルの出現確率の推定です。
90%信用区間（上下それぞれ５％）は、下側が3.28%上側が6.97%という結果です。
f:id:hippy-hikky:20190307003947p:plain

現在の期待値は5.05%です。

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

今のところの金のエンゼルの期待値は0.27%です。

計測結果

基礎集計

エンゼル出現確率の予測

広告

はじめに

資料

セミナー概要

正味重量のモデル

最尤推定とは

最尤推定の計算

ここまでの参考資料等

ベイズ推定

ベイズ推定とは

ベイズ推定で正味重量の推定

ベイズ推定で正味重量の推定結果

ここまでの参考資料等

階層ベイズモデル

階層ベイズモデルとは

階層ベイズモデルの実装

ここまでの参考資料等

まとめ

参考文献

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告

はじめに

問題設定

アプローチ

データ

データの概要

購入者毎のエンゼル出現確率の推定

最尤推定量

独立に出現確率を推定

購入者の個人差をパラメータに加えて推定

階層モデルとして予測

まとめ

参考資料

計測結果

基礎集計

エンゼル出現確率の予測

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告