【読後メモ】データ解析のための統計モデリング入門（その3）

改めて「データ解析のための統計モデリング入門」を読みました(一部飛ばしたけど)。本記事は読後のメモです。
全3回の連載で、今回が第3回目です。GLMのベイズモデル化(第9章)と階層ベイズモデル(第10章)についてまとめます。
第2回、第1回も合わせてご覧いただけたらと思います

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (29件) を見る

【目次】

はじめに
第9章 GLMのベイズモデル化と事後分布の推定
第10章階層ベイズモデル -GLMMのベイズモデル化-
まとめ
終わりに
参考文献

はじめに

「データ解析のための統計モデリング入門」を読んだので、どのようなことが書かれているのかといったメモを残しておきます。

この書籍では、統計モデリングとして「一般化線形モデル」に着目し、その基礎と発展について書かれています。よくある統計学の入門書では、確率の基礎や確率分布の紹介が多くの紙面を割いていることがあります。しかし本書では、一般化線形モデルに焦点を絞ることで、確率・統計の基礎についてはあまり触れられてはいません。これだけ書くと難しいように感じるかもしれませんが、逆に焦点を絞ることで本書を通して一貫した話題になっており、一般化線形モデルについて理解しやすいのではないかと思います。また、例題をベースにしているのも理解が進みやすいポイントかなと思います。

本書籍の読後メモは全3回でまとめる予定です。

本記事は第3回です。第3回では、第9章 GLMのベイズモデル化と第10章階層ベイズモデルについてまとめていきます。

【トップに戻る】

第1回 : learning-with-machine.hatenablog.com

第2回： learning-with-machine.hatenablog.com

第3回 : 本記事

第9章 GLMのベイズモデル化と事後分布の推定

第8章ではMCMC アルゴリズムを軸に、ベイズモデリングという考えの触りを学習しました。第9章では、第3章で学習したGLM（一般化線形モデル）のベイズ化を扱っています。ベイズ化/ベイズモデリングについては、なぜこれが必要なのかなど説明が本書籍ではあまり触れられていないように感じます*1が、この辺りについては私が書いた以下の記事などをご覧いただけると何かわかるかもしれません。

chocolate-ball.hatenablog.com

ざっくりと私の考えを述べると次のようになります。第7章までで扱ってきた最尤法ではパラメータを一つの推定値として求めるのに対し、ベイズ推定では事後分布という形でパラメータを推定します。最尤法は一点で予測するためにデータが少ない場合に極端な値を推定することがあります（過適合）。これに対してベイズ推定では、データが少ない場合には曖昧さを分布の広がりで表現できるということに違いがあると考えています*2。

ということで、第7章までで扱ってきた最尤法によるパラメータ推定をベイズで行うとどのようになるのかを見ていきます。なお、ベイズでのパラメータ推定には第8章で扱ったMCMC アルゴリズムを使用します。第8章では、種子の生存個数を二項分布としてその生存確率という単一のパラメータ推定問題を例として MCMC アルゴリズムが解説されていましたが、この章では複数のパラメータを推定する問題を扱います。

例題の設定

第4章で扱った植物の種子数のモデルが本章で扱う例題とされています(ポアソン回帰)。

個体 $i$ の種子数 $y_i$ は、パラメータ $\lambda_i$ のポアソン分布に従うと仮定します。ポアソン分布のパラメータ $\lambda_i$ は、個体サイズ $x_i$ を特徴量として予測されると仮定します。これらを表現すると以下のようなモデルになります。

${ \displaystyle y_i \sim Poi(\lambda_i) \\ \lambda_i = \exp(\beta_1 + \beta_2 x_i) }$

ここで、線形予測子のパラメータ $\beta_1, \beta_2$ が予測対象のパラメータで、これらをベイズモデリングにより推定していきます。なお第3章では、Rのglmライブラリを使って最尤推定しました。

GLMのベイズモデル化

上記のモデルをベイズモデル化します。ベイズの定理より、事後分布は以下の式で得られます。

${ \displaystyle p(\beta_1, \beta_2 | D) \propto p(D | \beta_1, \beta_2) p(\beta_1, \beta_2) }$

右辺第１項は尤度関数であり、パラメータ $\beta_1, \beta_2$ がある値の時のデータDの同時確率です。第２項はパラメータ $\beta_1, \beta_2$ の事前分布であり、本書籍では、 $\beta_1, \beta_2$ は独立と仮定しています( $p(\beta_1, \beta_2)=p(\beta_1) p(\beta_2)$ )。

最尤法では尤度関数だけを扱ってきましたが、ベイズ推定では事前分布を設定します。この事前分布がパラメータの範囲を定めるため、正則化をしているとみなせます（PRML 第1章）。次節から、尤度関数と事前分布を具体的に設定し、パラメータを実際に推定してみます。

尤度関数

尤度関数は、パラメータ $\beta_1, \beta_2$ がある値の時のデータDの同時確率です。ここで、データDはi.i.d.（独立同分布）を仮定するので、尤度関数は以下の式となります。

${ \displaystyle L(\beta_1, \beta_2 | D) = \prod_i p(y_i | \beta_1, \beta_2, x_i) }$

ここで、データDは $D=\{y_1, y_2, \cdots, y_n\}$ であり、 $x_i$ は植物iの個体サイズです。

事前分布

事前分布については、第8章で述べている通り、データがない時のパラメータの分布を表します。分析者の信念やドメイン知識に基づいて事前分布を設定することもありますが、本書籍では「無情報事前分布」を使うとしています。

無情報事前分布ですが、今回の推定対象であるパラメータ $\beta_1, \beta_2$ は、値の範囲が $[-\inf, \inf$ ]であるため、一律に特定の値をとるようにしてしまうと正規化して1になりません（発散する）。そのため、このような定義域のパラメータの無情報事前分布として、分散が大きな平たい正規分布を利用することが多いです。

事後分布の推定

上記の通り事前分布に平べったい正規分布を利用するとして、具体的に事後分布を推定していきます。実際の計算は、8章で述べられているMCMC アルゴリズムを利用します。

本書籍では、R言語とWinBUGSを利用したコードが紹介されています。 WinBUGS（他にもStan）とは、モデルとデータを定義することでMCMCサンプリングを実施してくれるソフトウェアのことで、 Probabilistic Programming Language（確率プログラミング言語、PPL）と呼ばれています。 WinBUGSやStanは、R言語とのインターフェースが用意されており、R界隈ではよく利用されていると思います。

しかし、私はRはあまり使っていないのでPythonを利用したいと思います。 Pythonに於けるPPLは、 PyMC3、TFP(Tensorflow Probability)、pyroなどが有名なところです（PyStanというものもありますが、私は詳しくありません）。以下のコード例はpymc3を利用したコードになります。

モデルの定義はセル番号5番で行なっています。ポアソン回帰なので、モデルとしては以下の通りです。

${ \displaystyle y \sim Poi(\lambda_i) \\ \log(\lambda_i) = \beta_1 + \beta_2 x_i }$

$\lambda_i$ は対数リンク関数で線形予測しているのですが、コード上では左辺の対数を右辺にて指数の形で表現しています。

パラメータ $\beta_1, \beta_2$ の事前分布は、ともに $N(0, 10)$ という正規分布としています*3。

以上のモデルを図で表現すると以下のようになります。

f:id:hippy-hikky:20190409001951p:plain — 種子数予測のためのモデル（Kruschkeのダイアグラム）

パラメータの予測結果は、以下のように事後分布で予測されています。

f:id:hippy-hikky:20190408235613p:plain — MCMCサンプリングの結果。上段が $\beta_1$ 、下段が $\beta_2$ であり、左列がヒストグラム（KDE）、右列がサンプル系列。

この結果から、事後分布の平均として、（ $\beta_1=1.29, \beta_2=0.08$ ）という推定結果を得ることができました。これは、第3章の最尤法での結果とほぼ一致します。

次に、事後分布からサンプルしたパラメータを利用して、回帰直線を描いてみます。

f:id:hippy-hikky:20190409000036p:plain — 推定した回帰直線。青点はデータ。オレンジ線が事後分布からサンプルした回帰直線で、青線は事後分布の期待値。

このようにして、GLMをベイズ化してパラメータが予測できました。

【トップに戻る】

第10章階層ベイズモデル -GLMMのベイズモデル化-

第9章ではGLM（ポアソン回帰）をベイズモデル化しました。続いて、第10章ではGLMMをベイズ化します。 GLMM（一般化線形混合モデル）については、当シリーズの第2回で触れていますので詳しくはこちらをご参照いただけたらと思います。

なおGLMMの最尤推定については、RのglmmMLパッケージなどがあるようですが、ちょっと調べたところPythonにはぴったりなパッケージが無さそうです。そのため、Pythonを使ってGLMMを解くためには、ここで紹介するように階層ベイズモデルにしてMCMC アルゴリズムで推定するのが簡単ではないかなと思います。

例題の設定

第7章で扱った個体差のある種子の生存確率のモデリングが例題とされています。「各個体で8個の種子をサンプルし、その中で何個の種子が生存していたか」という問題のため、普通に考えると二項分布でモデリングするでしょう。しかし、個体差があるため、下図の通り過分散なデータになっています。

f:id:hippy-hikky:20190409215822p:plain — 例題データのヒストグラム。各個体8個の種子をサンプルし、生存していた種子の個数についてのヒストグラム。

第7章では、過分散が個体差によるものとして、GLMMを使ってモデリングしてきました。これを階層ベイズモデルとしてモデリングしてみます。

階層ベイズモデル

8,9章で扱ったベイズモデルでは、求めたいパラメータに事前分布を設定しました。階層ベイズモデルとは、事前分布を構成するパラメータもデータから推定するためにさらに事前分布（階層事前分布と呼ぶ）を設定します。この階層モデル（階層事前分布）の意味を例題に従ってみていきます。

まず本データである個体 $i$ の生存種子数 $y_i$ は第7章と同様に以下の二項分布に従うと考えます。

${ \displaystyle y_i \sim \binom{8}{y_i}q^{y_i}_{i}(1-q_i)^{8-y_i} }$

ここで、二項分布のパラメータである生存確率 $q_i$ ですが、全ての個体で同じ生存確率という単純な二項分布で表現すると上記の通り個体差があるためモデルの精度は悪化します（過分散）。そこで、 $q_i$ は以下で示すように、個体差 $r_i$ を持つとしてみます（第7章）。

${ \displaystyle logit(q_i) = \beta+r_i }$

$\beta$ は全個体に共通のパラメータです。個体差が少ない均質な植物群であれば、個体差rは小さく、単純な二項分布モデルと同じものとして表現できます。第7章では、個体差 $r_i$ を積分消去して期待値を計算していましたが、階層ベイズモデリングではこの $r_i$ にも事前分布を設定して、推定します。

まず、 $r_i$ は $N(0, s^{2})$ の正規分布に従うものとします（第7章）。ここに出てくる $s$ は、個体毎のばらつきの大きさを規定するパラメータです。 $s$ が小さければ、植物群は全体的に均質（ $r_i$ は0に近いものばかり）であると見れるし、 $s$ が大きければ、植物群のばらつきが大きいと言えます。このように、推定対象のパラメータ $r_i$ の全体の動きを規定するための（上位階層の）パラメータを仮定して、一緒に推定するモデルが階層ベイズモデルです。

なお、個体差 $r_i$ の事前分布のように、事前分布の事前分布を階層事前分布と呼び、階層事前分布のパラメータ（ばらつきの大きさ $s$ ）をハイパーパラメータと呼びます。

階層ベイズモデルの推定

では、上記のモデルを使って、具体的にパラメータの推定をしてみます。先の9章での実験の通り、PyMC3を利用して計算してみます。 jupyter notebookを添付しますので、ご参照ください。

モデルの定義はセル番号5で行なっています。先に示した通り、種子数は二項分布に従うとし、そのパラメータである生存確率 $p_i$ を共通パラメータ $\beta$ と個体差 $r_i$ で表現しています。

$\beta$ の事前分布は平たい正規分布としています。個体差 $r_i$ の事前分布は $N(0, s^{2})$ の正規分布です。そして、 $s$ の事前分布として一様分布を設定しています。

以上のモデルを図で表現すると以下のようになります。

f:id:hippy-hikky:20190410002523p:plain — 生存種子数のモデル（Kruschkeのダイアグラム）

個体差 $r_i$ 全体を統括するハイパーパラメータを、個体に関係なく共通のパラメータ $s$ として定義します。これは先に記したように、植物群の個体差がどの程度あるのかを規定するパラメータです。

各パラメータの事後分布は以下の通りに推定されました。

f:id:hippy-hikky:20190410002721p:plain — MCMCサンプリングの結果。上段が個体差の大きさを示すパラメータ $s$ 、中段が各個体の $r_i$ 、下段が共通のパラメータ $\beta$ 。左列がヒストグラム（KDE）、右列がサンプル系列。

MCMCサンプリングの結果。上段が個体差の大きさを示すパラメータ $s$ 、中段が各個体の $r_i$ 、下段が共通のパラメータ $\beta$ 。左列がヒストグラム（KDE）、右列がサンプル系列。

$\beta$ と $s$ の事後分布を見てみると、書籍内の事後分布と概ね同じような値であることがわかります。この結果（事後分布）を利用して、生存種子数のシミュレーションをした結果を以下に示します。

f:id:hippy-hikky:20190410003142p:plain — 生存種子数のシミュレーション結果。青線が元のデータ。オレンジが推定値。

推定した事後分布を利用して100個体それぞれから10セットサンプルし、平均値をプロットしています。推定値がデータとよく合っているであろうことがわかります。

ベイズモデルで使う事前分布

本書籍ではここまで、「事前分布」についてはあまり深く説明されてきませんでした。本記事では、事前分布とはデータが無い時のパラメータに対する信念を表現するものだと書きましたが、この信念が大きく間違っていると推定に悪影響が現れます。極端な例では、事前分布で確率を0とした値の範囲には事後分布が現れません*4。このような極端な例ではなくても、事前分布がずれていると、合理的な事後分布を推定するためにたくさんのデータが必要になってしまいます。

本書籍では、このように事前分布の選択はベイズモデリングにおいて重要であるため、主観的な事前分布は使う必要がないと述べてられています。基本的には無情報事前分布を使えば良いと書かれています。しかし、上記の種子の生存確率の個体差のような問題の場合、個体差 $r_i$ の事前分布に無情報事前分布を利用するということは、各個体のそれぞれを独立に推定していることになるため、統計モデリングとしては望ましくありません。そのため、このような場合には階層事前分布を設定することモデルとしては好ましいでしょう（と書籍内に記載されています）。

階層事前分布を使用するべきなのか無情報事前分布を使用するべきなのかについては、本書籍で「パラメータがデータのどの範囲を説明するのかに依存する」と主張されています。同じようなパラメータが少数で大域的に影響するパラメータ（上記の $\beta$ のようなパラメータ）には無情報事前分布を使い、同じようなパラメータが多数で局所的に影響するパラメータ（上記の個体差 $r_i$ ）には階層事前分布を利用するべきとされています。局所的なパラメータは、独立に推定してしまうと、データに過適合してしまう恐れがありますから、階層事前分布で正則化することは尤もなやり方でしょう。

なお私は、主観的な事前分布も場合によっては使うべきだと考えています。データが少ない場合で、且つ、問題に対しての深い知識/常識的な知識がある場合には、積極的にパラメータの範囲を拘束するべきと考えています。データの範囲を拘束することで、データが少なくてもリーズナブルな事後分布を得られる場合があるためです*5。

【トップに戻る】