PyMC3 を使ったベイズ一般化線形モデル回帰分析メモ

久々にIPAソフトウェア開発データ白書のデータを分析してみようと思い、ベイズ推定を使ってみたのだけど、データの分散が大きすぎるので分位点回帰みたいな方法使わないとだめだということがわかったのでメモだけ。

インストール

環境は Windows 10 の WSL2 Ubuntu 20.04 LTS。WSL2 は本当に便利。

sudo apt install python3 python3-pip python3-pandas jupyter-notebook
pip3 install pymc3

WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions. という Waning が出る場合は、MKL という NumPy を高速化する Intel 謹製ライブラリをインストールすると消える。

pip3 install mkl

ライブラリとデータの読み込み

ここは例によって例の通り。以下、y、x1、x2 の三列がある前提で書いている。

import pandas as pd
import matplotlib.pyplot as plt
import pymc3 as pm

indata = pd.read_csv('./indata.csv') # y, x1, x2

ベイズGLM

GLM用の便利関数が用意されているし、formula も使えるので一見簡単に使える。 これが実際には、いろいろと初期設定を変える必要があり結構詰まる。

with pm.Model() as model:
    pm.glm.GLM.from_formula('y ~ x1 + x2', indata)
    trace = pm.sample(draws=5000, chains=4, cores=1, start=pm.find_MAP(), step=pm.NUTS())

formula は indata の中の変数名を使える。y ~ x なら単回帰、y ~ x1 + x2 なら重回帰、切片はデフォルトであり扱いなので不要な場合は y ~ x1 + x2 -1 のように -1 をつけると切片=0になる。

draws はサンプルの生成回数。chains は生成するサンプル列の数。ベイズ推定の場合、MCMCを使うためいつも異なる結果を返すので、複数回実行して安定した結果を得られているか確認することができる。cores は並列実行数。

start は、初期の事前分布を設定するもので、find_MAP は MAP推定値を設定するものらしい。正しい使い方なのかはよくわからない。

step は、サンプルを生成するサンプラー。2値の場合は BinaryMetropolis、離散値の場合は Metropolis、連続値の場合は NUTS を使うと良いらしい。

from_formula メソッドの引数に family でリンク関数を priors を付けることで説明変数の事前分布を変えることができる。family に指定できるリンク関数は、Normal（正規分布）、StudentT（スチューデントのt分布）、Binomial（2項分布）、Poisson（ポアソン分布）, NegativeBinomial （負の2項分布）の5つ。priors で使える分布はここで説明されているものが利用できる。後者はグラフとともに説明されているので大変わかりやすい。

    pm.glm.GLM.from_formula('y ~ x1 + x2', indata, 
         family = pm.glm.families.Normal(),
         priors = {
            'x1': pm.Normal.dist(mu=0., sigma=1.),
            'x2': pm.Normal.dist(mu=0., sigma=1.),
            'Intercept': pm.Normal.dist(mu=0., sigma=1.), // 切片は Intercept
        }
    )

trace にはサンプリングの結果が入っている。この結果は PyMC3 付属のプロット系関数でグラフ化して確認できる

pm.traceplot(trace[4000:]);
pm.summary(trace[4000:])
pm.plot_posterior(trace[4000:]);

trace[4000:] と書いているのは、5000回サンプルを生成したうち、安定したであろう 4000回以降を使うという意味。

traceplot でデータが収束しているかを判断できる。左側のグラフは同じ形で重なっている、右側は特定の値周辺に散らばっていれば、収束しているとみなしていい。summary の r_hat が 1.1 以下という判断でも良いらしい。

summary は各値の詳細、plot_posterior は summary 内容を可視化したもの。