2019-04-03

Chainer Meetup #09 でLT：チョコボール数の傾向を分析

チョコボール chainer 分析 Presentation

【概要】

3/30に開催のChainer Meetupに参加してきました
ただ参加するだけではつまらないので、LTしてきました
本記事ではLTの概要と資料をまとめます

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

【目次】

はじめに
資料
LT概要
データと分析アプローチ
- データ収集効率化の取り組み
- 実験に利用したデータの概要
モデル
結果
終わりに
参考文献
広告

はじめに

chainer.connpass.com

今回で2回目の参加となりましたChainer Meetupですが、今回もとても興味深い話が聞けて良かったです。 Chainerの中の人の開発情報が聴けるというのはもちろん良かったのですが、 nVidiaの中の人の講演が今回もとても興味深かったです。前回もパフォーマンスチューニングについてとても興味深く聞いていましたが、今回は自分でも取り組めそうな話で希望が持てました。資料はconnpassのイベントページに公開されているので、興味ある方はご覧になられたらと思います。

本記事ですが、前回に引き続き今回もLTをさせてもらったので、どのような話をしたのか概要をまとめておくものです。

【トップに戻る】

資料

当日使用した資料は以下です。

20190330 chainer meetup chocoball from Yoichi Tokita

【トップに戻る】

LT概要

今回のLTでは、チョコボールの個数を分析した話をしました。

チョコボール（ピーナッツ味）は、内容量が28gと定められています*1。内容量（重量）が定められているのですが、ピーナッツは自然のものであるため大きさや重さを厳密に一定にすることはできません*2。そのため、一箱に入っているチョコボールの個数は箱ごとに異なります。

一方、せっかくチョコボールを買うなら美味しいものを食べたいのですが、小さかったり形がいびつなチョコボールは味もいまいちな印象です(個人的な感想であり品質を保証するものではありません)。

そこで、チョコボールの粒が大きく揃っているものを購入するにはどのようにすれば良いのかについて、個数の傾向から考察してみたという話です。

【トップに戻る】

データと分析アプローチ

データ収集効率化の取り組み

（詳しくは上記の資料を参照ください）

個数の傾向を分析するには、各商品に個数がいくつ入っていたのかを計測し、データを利用できる形に整理することが必要です。しかし、一つ一つ数えるのはなかなかの手間です*3。

そこで、Faster R-CNNを使ってチョコボールを検出するWebサービスを開発しました。これで、写真を撮るだけで個数が自動で計測できます！（ここまでは、前回のLTでお話ししたところです）

www.chocolate-ball.net

しかしながら実際のところ、このシステムは計測の場面ではまだ実運用できておらず、手で数えています*4。 ChainerのイベントでChainerについて一切触れないのはマズかろうということで(大人の事情)。

実験に利用したデータの概要

ということで手で数えて収集したデータですが、2017年11月から2019年3月まででピーナッツ味373箱分のデータを利用します。

f:id:hippy-hikky:20190401222052p:plain — チョコボール（ピーナッツ味）の個数のヒストグラム

このデータを賞味期限が同じチョコボールをひとまとめにして、個数の分布（箱ひげ図）を確認してみてみます。なお、チョコボールには製造年月日の記載が無く、賞味期限が月単位で表示されています。そのため、賞味期限が同じであれば、同じ時期に製造されたものであると仮定して、データを分析します。

f:id:hippy-hikky:20190401222325p:plain — 賞味期限ごとの個数分布。横軸はパッケージに記載の賞味期限。縦軸はチョコボールの個数。

上記の図からわかるように、2018年の8~10月の個数が少ないように見えますね。また、2019年の8~10月もデータ数が少し少ないですが、同様に個数が少ないように見える気がします。

そのため今回は、賞味期限が8~10月の個数が他の期間と比べて有意に数が少ないのかを確認してみます。この期間の数が少ないことが統計的に確認できたならば、 賞味期限が8~10月のチョコボールを狙って購入することで、粒の大きいチョコボールを手に入れられる可能性が高まる と期待できます。

他にも分析する軸は考えられますが、今回は時期による個数の変化にだけ注目します。

【トップに戻る】

モデル

個数の傾向を分析するに当たって、下記の統計モデルを設定しました*5。

f:id:hippy-hikky:20190401223413p:plain — チョコボールの個数に関しての統計モデル（階層ベイズモデル）。

一箱に入っているチョコボールの個数yは、ポアソン分布 $Poi(\lambda)$ に従う確率変数と考えます。ポアソン分布のパラメータ $\lambda_i$ は各月 $i$ の平均個数を示すパラメータであり、ガンマ分布 $Gamma(\alpha, \beta)$ を事前分布とします。ガンマ分布のパラメータ $(\alpha, \beta)$ は任意に設定するのではなく、それぞれ階層事前分布を設定してデータから推定させることにします。ガンマ分布のパラメータ $(\alpha, \beta)$ は、平均 $\mu$ と標準偏差 $\sigma$ を示すパラメータを使って以下のように算出されます。

$\displaystyle{ \alpha = \mu^{2}/\sigma^{2} \\\ \beta = \mu / \sigma^{2} }$

今回は、8~10月の期間のデータから推定した事前分布( $Gamma(\alpha_0, \beta_0)$ )と、その他の期間のデータから推定した事前分布( $Gamma(\alpha_1, \beta_1)$ )を比較し、統計的に差があるのかを確認します。

【トップに戻る】

結果

上記のモデルを使い、パラメータを推定しました。
パラメータは、MCMC アルゴリズム（NUTS）を利用して数値的に推定しました。

推定した二つのガンマ分布は以下のようなものです。

f:id:hippy-hikky:20190401223516p:plain — 事前分布の推定結果。上が8~10月の期間の推定事前分布。下がその他の期間の推定事前分布。8~10月の期間が平均で約1.4個少ない。

95%信用区間で分布の平均に差があるのかを確認するために、平均のパラメータの事後分布の差の分布を以下に示します。

f:id:hippy-hikky:20190401223705p:plain — 推定された事前分布のパラメータの差の事後分布。左が平均の差で、90%信用区間が0を超えており、差があることを示している。右は分散の差で、こちらはモードが0の付近であり、差がないことがわかる。

ということで、二つの分布に差がないことを示す0の位置が90%信用区間から外れているため、二つの分布の平均には差があるであろうということがわかりました。

その差は、約1.3個。

以上の結果、賞味期限が8~10月のチョコボールは個数が少ないということがわかりました。資料p.22にポアソン分布の各月のパラメータ $\lambda_i$ の事後分布を載せていますが、 2018年だけでなく、2019年の8~10月も同様に個数は少ないと言えそうです。そのため、賞味期限が8~10月のチョコボールを狙って買えば、粒の大きいチョコボールを手に入れられそうです！

【トップに戻る】

終わりに

ということでチョコボールの個数の傾向を分析し、特定の月の個数に差がありそうだということまではわかりました。

しかし、今回は8~10月を予め切り分けてパラメータを推定したので、 8~10月だけが少ないのかはわからないです。データから変化点を合わせて推定させることが必要と考えています*6。
（実は、変化点を推定する実験を最初にやってみたのですが、推定結果が安定せずでして、プログラム上に問題があるのか、モデルがまずいのか検討中の状態です。。。）

また、チョコボールの大きさが大きくて粒が揃っているというのは、個数だけで特徴付けられるのかについて、はっきりとした確証はありません。そのため、大きさのばらつきについても今後調査していきたいと考えています。

【トップに戻る】

参考文献

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

【トップに戻る】

第143回チョコボール計測（チョコバナナ）

計測チョコボール

本日の計測結果を報告します。今回からチョコバナナ味の計測が続きます。

チョコバナナ味の重量傾向を調査したかったので、在庫がなくなる前にAmazonで購入しておきました。

今日は2箱計測します。

計測結果

date	best_before	weight	box_weight	number	factory	shop	angel	net_weight	mean_weight
2019-03-27	2019-11-01	27.829	4.859	15	小山工場	amazon	なし	22.970	1.531
2019-03-27	2019-11-01	27.665	4.854	15	小山工場	amazon	なし	22.811	1.521

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

基礎集計

チョコバナナ味の集計結果です。

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

	最小	中央値	最大値	平均
正味重量	22.107	22.599	24.338	22.759
個数	14.000	15.000	16.000	14.800

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

エンゼル出現確率の予測

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

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

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

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

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

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

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

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

計測チョコボール

本日の計測結果を報告します。今回も通常版のピーナッツを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-26	2019-08-01	34.040	4.732	16	小山工場	スーパー（さいたま市）	なし	29.308	1.832
2019-03-26	2019-08-01	34.223	4.707	15	小山工場	スーパー（さいたま市）	なし	29.516	1.968
2019-03-26	2019-08-01	35.153	4.760	15	小山工場	スーパー（さいたま市）	なし	30.393	2.026
2019-03-26	2019-08-01	34.214	4.742	16	小山工場	スーパー（さいたま市）	なし	29.472	1.842

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

基礎集計

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

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

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

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

エンゼル出現確率の予測

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

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

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

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

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

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

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

第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%です。

はじめに

資料

LT概要

データと分析アプローチ

データ収集効率化の取り組み

実験に利用したデータの概要

モデル

結果

終わりに

参考文献

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告

はじめに

資料

セミナー概要

正味重量のモデル

最尤推定とは

最尤推定の計算

ここまでの参考資料等

ベイズ推定

ベイズ推定とは

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

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

ここまでの参考資料等

階層ベイズモデル

階層ベイズモデルとは

階層ベイズモデルの実装

ここまでの参考資料等

まとめ

参考文献

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告

計測結果

基礎集計

エンゼル出現確率の予測

広告