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

【概要】

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を利用した実装例もあるが、実装部分は無視しても理解できる。