PyStanで実装するベイズ推定（単回帰・ロジスティック回帰）


# Stan, PyStanとは

**Stan**とは、C++をベースに実装された確率的プログラミング言語です。NUTSアルゴリズム (HMCを発展させたもの) を用いて、様々なベイズ推定を行うことができます。特徴として以下のようなものが挙げられます。

* 統計モデルの記述が簡単
* HMCなのでサンプリングが高速
* さまざまな確率分布を利用可能

**PyStan**はStanをPythonから扱うためのインターフェースを提供するパッケージです。Stanの文法に従ってモデルを記述し、コード実行時にモデルをコンパイルします。

* [PyStan: The Python Interface to Stan](https://pystan.readthedocs.io/en/latest/)

# PySTanのインストール

pipでインストール可能です。トレースプロットの描画のために `arviz` もインストールしておくと良いです。

```
$ pip install pystan
$ pip install arviz
```

# PyStanによるベイズ推定

それでは実装していきます。まず必要なライブラリをインポートします。

```python:In
import numpy as np
import pandas as pd
from scipy.stats import mstats

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

import pystan
import arviz
```

## 1. 単回帰

ここではAuto MPGデータセットを使用します。[UCI Machine Learning Repository](https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/) からダウンロードできます。
データをpandas DataFrameに格納します。

```python:In
column_names = ['MPG','Cylinders','Displacement','Horsepower','Weight',
                'Acceleration','Model Year','Origin']
df = pd.read_csv('auto-mpg.data', names=column_names, na_values = "?", 
                 comment='\t', sep=" ", skipinitialspace=True)
```

燃費 (`MPG` : Mile per Gallon) は車両の重量 (`Weight`) に反比例するのではないか、という予想の下でこれを確かめてみます。

```python:In
df['Weight_inv'] = df['Weight'].apply(lambda x: 1/x)
plt.scatter(df['Weight_inv'], df['MPG'])
```

![mpg_check.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/304b184e-1381-929b-6c8e-2e51a198b03d.png)

確かに重量の逆数は燃費と正の相関にありそうです。横軸の値によってばらつきが異なっているようですが、ここではひとまず以下のような統計モデルを仮定します。

$$MPG \sim Normal(a × Weight^{-1} +b, \sigma)$$

すなわち燃費は、重量の逆数の線形関数を平均とし、$\sigma$を分散とする正規分布に従うと仮定します。ここで出てきた$a, b, \sigma$が推定の対象となるパラメータになります。この統計モデルをStanの文法に従って記述すると以下のようになります。

```python:In
stan_model = """
data {
    int N;
    real X[N];
    real Y[N];
}
parameters {
    real a; 
    real b;
    real<lower=0> sigma;
}
model {
    for (n in 1:N) {
        Y[n] ~ normal(a * X[n] + b, sigma);
    }
}
"""
```

Stanのコードは以下のブロックが必要になります。

* **data** ... 観測データを記述
* **parameters** ... 推定対象のパラメータを記述
* **model** ... 仮定した統計モデルを記述

分散 `sigma` は非負の実数値をとるので、データ型を `real<lower=0>` とします。
モデルの記述が完了したら、これをコンパイルします。

```python:In
sm = pystan.StanModel(model_code=stan_model)
```

観測データを辞書型で与え、MCMCサンプリングを行います。データ辞書のkeyは、Stanコードのdataブロックで定義した変数名に対応します。
`pystan.StanModel.sampling` の引数として、`iter` と `warmup` を与えています。`iter` はモンテカルロステップ数です。`warmup` ですが、定常状態に至るまでにしばらくかかるため、初期段階のサンプリングの結果は捨てた方が良いとされます。最初の何ステップを捨てるかを `warmup` で指定しています。

```python:In
stan_data = {
    'N': df.shape[0], 
    'X': df['Weight_inv'], 
    'Y': df['MPG']
}
fit = sm.sampling(data=stan_data, iter=2000, warmup=500, chains=3, seed=123)
print(fit)
```

```:Out
Inference for Stan model: anon_model_cdf2dd4817e02d2c4cf1f2106bb8e4c0.
3 chains, each with iter=2000; warmup=500; thin=1; 
post-warmup draws per chain=1500, total post-warmup draws=4500.

        mean se_mean     sd   2.5%    25%    50%     75%  97.5%  n_eff   Rhat
a      6.6e4   54.18 2126.0  6.2e4  6.5e4  6.6e4   6.8e4  7.0e4   1540    1.0
b      -0.57    0.02   0.81  -2.12  -1.13  -0.56 -8.4e-3   1.01   1576    1.0
sigma   4.25  3.5e-3   0.16   3.94   4.14   4.24    4.35   4.56   1962    1.0
lp__  -772.7    0.03   1.24 -775.9 -773.3 -772.4  -771.8 -771.3   1571    1.0

Samples were drawn using NUTS at Sat Dec 19 07:48:35 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
```

各推定値のRhatが1.0となっているため、定常状態に落ち着いていると解釈できます。
また、各推定値の平均値、50%信頼区間、95%信頼区間を確認できます。

各推定値のサンプリング結果は以下のようにアクセスできます。

```python:In
fit.extract('a')['a']
# array([66799.22277616, 62788.84373091, 63577.82392938, ...,
#        63620.47069878, 63571.6229788 , 67382.64630103])
```

サンプリング結果の分布とトレースプロットを次のように表示できます。

```python:In
fig = arviz.plot_trace(fit)
```

![スクリーンショット 2020-12-19 8.44.53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/2e840e96-1947-085e-5059-d88557598d8b.png)

それではa, bの推定値を用いて予測直線を信頼区間とともに描いてみましょう。

```python:In
a, b = 6.6e4, -0.57
x = np.linspace(0.0002, 0.0006, 100)
y = a * x + b
plt.plot(x, y, color='red')
plt.scatter(df['Weight_inv'], df['MPG'])
```

![mpg_result.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/bdbcecfa-778a-2489-5d6d-628c8e32fc3d.png)

概ね良い推定が出来ていることが確認できます。

## 2. ロジスティック回帰

続いてロジスティック回帰です。今度は予測回帰曲線を信頼区間とともに表示してみましょう。
scikit-learnの乳がんデータセット (breast_cancer) を使用します。細胞の半径 (`mean radius`) を説明変数、悪性かどうか (`malignant`) を目的変数とします。

```python:In
from sklearn.datasets import load_breast_cancer

bc = load_breast_cancer()
df = pd.DataFrame(bc.data, columns=bc.feature_names)
df['benign'] = bc.target
df['malignant'] = df['benign'].apply(lambda x: 1 if x == 0 else 0)
```

散布図を確認します。

```python:In
plt.scatter(df['mean radius'], df['malignant'])
plt.xlabel('Mean radius')
plt.ylabel('Malignant')
plt.xlim([5, 30])
plt.yticks([0, 1])
```

![bc_check.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/88b01e50-edb1-1724-0036-a2f9ac233feb.png)

細胞の半径 (`mean radius`) が大きいほど、悪性である確率が大きい傾向が読み取れます。以下のようなロジスティック回帰モデルを仮定します。

$$Malignant \sim \exp(a \times MeanRadius + b))$$

単回帰同様、Stanコードを記述し、コンパイルします。
モデルの記述において、ロジスティック回帰では `bernoulli_logit` を使用します。

```python:In
stan_model = """
data {
    int N;
    real X[N];
    int<lower=0, upper=1> Y[N];
}
parameters {
    real a;
    real b;
}
model {
    for (n in 1:N) {
        Y[n] ~ bernoulli_logit(a*X[n] + b);
    }
}
"""
```

観測データを与え、MCMCサンプリングを行います。
また、トレースプロットとともに出力結果を表示します。

```python:In
stan_data = {
    'N': df.shape[0], 
    'X': df['mean radius'], 
    'Y': df['malignant']
}
fit = sm.sampling(data=stan_data, iter=3000, warmup=500, chains=3, seed=123)
print(fit)

fig = arviz.plot_trace(fit)
```

```:Out
Inference for Stan model: anon_model_961ca63f5fca60aa2e92ea6e3198f2fb.
3 chains, each with iter=3000; warmup=500; thin=1; 
post-warmup draws per chain=2500, total post-warmup draws=7500.

       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
a      1.05  2.8e-3   0.09   0.88   0.98   1.04   1.11   1.24   1166    1.0
b    -15.46    0.04   1.35 -18.23 -16.32 -15.41  -14.5 -13.02   1160    1.0
lp__ -166.0    0.02   0.98 -168.6 -166.4 -165.7 -165.2 -165.0   1896    1.0

Samples were drawn using NUTS at Sat Dec 19 13:40:51 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
```

![スクリーンショット 2020-12-19 14.04.32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/b772b514-ea35-30af-25a7-06615495b778.png)

RHat=1.0ということで、よく収束していることが分かります。
a, bの推定値は回帰曲線の表示に利用するので変数に格納しておきます。

```python:In
a, b = 1.05, -15.46
```

続いて、サンプリング結果を用いて信頼区間を求めていきます。
前述の通り、`extract()` メソッドでサンプリング結果を抽出します。予測区間を [0, 30] としてサンプリング回数 (2500×3=7500回) 分のロジスティック回帰曲線を求め、pandas DataFrameに格納します。計算結果をもとに、`scipy.stats.mstats.mquantiles` を用いて50%信頼区間、99%信頼区間を求めます。

```python:In
ms_a = fit.extract('a')['a']
ms_b = fit.extract('b')['b']

x = np.linspace(0, 30, 100)
logistic_func = lambda x: 1.0 / (1.0+np.exp(-x))
df_lines = pd.DataFrame([])
for i in range(x.shape[0]):
    df_lines[i] = logistic_func(ms_a*x[i] + ms_b)

low_y50, high_y50 = mstats.mquantiles(df_lines, [0.25, 0.75], axis=0)
low_y99, high_y99 = mstats.mquantiles(df_lines, [0.005, 0.995], axis=0)
```

それでは信頼区間とともに、観測データおよび回帰曲線を描画してみましょう。

```python:In
plt.fill_between(x, low_y50, high_y50, alpha=0.6, color='darkgray')
plt.fill_between(x, low_y99, high_y99, alpha=0.3, color='gray')
plt.plot(x, logistic_func(a*x+b))
plt.scatter(df['mean radius'], df['malignant'])
plt.xlabel('Mean radius')
plt.ylabel('Malignant')
plt.xlim([5, 30])
plt.yticks([0, 1])
```

![bc_display.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/bf3c9814-a0cb-2914-23cc-e08ad52c5b52.png)

# まとめ

以上、PyStanの紹介とPyStanによる単回帰とロジスティック回帰の実装でした！
良いデータセットが見つかれば階層ベイズなんかも記事に出来たらと思います。