pymcを使ってベイズ推定による線形回帰


## はじめに
　[Pythonで動かして学ぶ! あたらしいベイズ統計の教科書](https://www.amazon.co.jp/dp/4798168645)を読んだ事をgoogle colaboratryでやろうと思ったら[こちらの記事にあるように](https://discourse.pymc.io/t/google-colab-update-pymc3-is-being-replaced-by-pymc-v4-faq/10235)、colabではpymc3ではなくpymcを使ってくれとのことなので、pymcでやってみました。

## ベイズ推定とは
　ベイズ推定は、統計的な手法の1つであり、事前に与えられた確率分布（事前分布）とデータを元に、事後分布を計算し、その分布を利用してパラメータの推定や予測を行う手法です。

　ベイズ推定では、あるパラメータ（例えば、母平均や母分散など）が持つ確率分布を事前分布として仮定します。この事前分布は、事前にパラメータについて知っている情報や仮定に基づいて設定されます。そして、データが得られた後、データと事前分布を組み合わせて事後分布を計算し、それによりパラメータの推定や予測が可能になります。

## ベイズ推定に使用するデータの取得と前処理

　今回も、[私の以前書いた記事](https://qiita.com/miso_taku/items/552f75877014ef34d233)で使用したScikit-learnのデータセットのボストン市の住宅価格を使おうと思ったら、こいつもScikit-learnの最新バージョンじゃdeprecatedになってるじゃん...、仕方ないので、同じScikit-learnのデータセットにある糖尿病患者の診療データセットを使用したいと思います。糖尿病患者の診療データセットは、糖尿病患者の検査数値と1年後の進行状況をデータセットにしたものです。

 まずは、データ取得と前処理に使うライブラリのインポート

```python
import numpy as np
import pandas as pd
from sklearn.datasets import load_diabetes
```

　そしてデータの読み込み

```pyhton
diabetes = load_diabetes()
diabetes_data = pd.DataFrame(diabetes.data, columns = diabetes.feature_names).assign(target=np.array(diabetes.target))
```

　せっかくなので、データの中身と概要も見ておきましょう。

```Python
diabetes_data.head()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/350321/89e4c1f8-597f-68f5-1146-e453f4265fda.png)

```Python
diabetes_data.describe()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/350321/b945ed19-8707-1ddf-57b5-775cbbe24b37.png)

seabornのpairplotを使ってデータの関係性も見ておきましょう。

```python
import seaborn as sns

sns.pairplot(diabetes_data)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/350321/05287b86-1ea6-129a-fca5-f3d121081c39.png)

　とりあえずこの中で、s1(tc、総血清コレステロール)とs2(ldl、低密度リポタンパク質)が高い相関関係を持ってそうなので、それを使って単回帰分析を行います。

```python
X = diabetes_data[['s1']].values
y = diabetes_data[['s2']].values
```

 Scikit-learnを使った単回帰分析

```python
from sklearn.linear_model import LinearRegression

reg = LinearRegression()
reg.fit(X, y)
y_pred = reg.predict(X)
```

回帰直線の描画、単回帰分析なのでy=α＋βXで表せます。

```python
import matplotlib.pyplot as plt

plt.scatter(X, y)
plt.plot(X, reg.predict(X), color='black')
plt.xlabel('s1')
plt.ylabel('s2')
plt.grid()
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/350321/aa8cadab-b895-8c97-9251-e17f33442928.png)

 Scikit-learnによる単回帰分析で求めた回帰直線の切片(alpha_0)と傾き(beta_0)を表示

```python
alpha_0 = reg.intercept_[0]
beta_0 = reg.coef_[0][0]
print("alpha_0:", alpha_0)
print("beta_0:", beta_0)
```

出力
```
alpha_0: 5.1571733982230335e-17
beta_0: 0.8966629578104912
```

## pymcを使ったベイズ推定による線形回帰

　それでは、ここからが本題のpymcを使ったベイズ推定による線形回帰になります。

　pymcをインポートしてから、取得したデータを、切片(α)、傾き(β)の事前分布が平均0、標準偏差10の正規分布、観測値の標準偏差が半正規分布の線形回帰モデルを定義します。

```python
import pymc as pm

with pm.Model() as model:

    alpha = pm.Normal("alpha", mu=0, sigma=10) #切片
    beta = pm.Normal("beta", mu=0, sigma=10) #傾き
    sigma = pm.HalfNormal("sigma", sigma=1) #データのばらつき

    mu = alpha + beta * X

    Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=y)
```

　sample関数を用いて事後分布のサンプリング

```python
with model:
    trace = pm.sample()
```

　plot_trace関数を用いて結果を可視化

```python
var_names = ["alpha", "beta", "sigma"]

pm.plot_trace(trace, var_names=var_names, backend_kwargs={"constrained_layout":True})
```

出力
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/350321/b68f0f4a-52ba-aad8-1a7f-d663bfe0ed85.png)

　summary関数で要約統計量を表示。結果からScikit-learnによる単回帰分析で求めた回帰直線の切片と傾きとほぼ同じであることが確認できました。

```python
pm.summary(trace, var_names=var_names)
```

 出力
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/350321/af96345d-8262-b54f-a5d6-abb86b7dd6df.png)


　plot_posterior関数で推論された事後分布と信頼区間を表示。

```python
pm.plot_posterior(trace, point_estimate='mode')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/350321/b0dfbb4e-1c14-2877-dbbd-349d9b4b6f17.png)


## おわりに

　本当は以下みたいな感じの線形回帰のグラフに信頼区間も表示したかったのですが、pymc3からpymcに変えた影響もあってか、上手くできませんでした。機会があればまた挑戦したいと思います。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/350321/f8b7f8a0-0bb8-733a-f53f-28ec1e0dbaa4.png)