PyMC3を使ったPythonベイズ推定超入門（その3）ベイズ型線形回帰モデルを作ってみよう！

前回と前々回にベイズ推定について簡単に説明しました。

以下、前回の記事です。

以下、前々回の記事です。

通常の線形回帰モデルは、切片や係数、予測の値は1つです。

一方，ベイズ線形回帰モデルで得られるのは分布です。正確には、分布から得られるサンプルです。

今回は、PyMC3を使ってベイズ線形回帰モデルを構築する方法を紹介します。

通常の線形回帰モデルと比較します。

Contents

通常の線形回帰とベイズ線形回帰
データセット
準備
通常の線形回帰モデル
ベイズ線形回帰モデル
　パラメータの事後分布
　予測モデル
まとめ

通常の線形回帰とベイズ線形回帰

簡単に数式で考えていきます。

通常の線形回帰モデルは、以下のように表現されます。yが目的変数で、xが説明変数です。

$\displaystyle y = a + b x + \epsilon , \epsilon \sim \mathcal{N}(0,\sigma^2)$

これは、 $y$ が平均 $ax+b$ と標準偏差 $\sigma$ で正規分布していることを意味します。

$\displaystyle y \sim \mathcal{N}(a + b x,\sigma^2)$

この $a$ と $b$ は、ある固定した値です。

一方、ベイズ統計学では $a$ と $b$ を固定したものとして扱うのではなく、変化するものとして扱います。

例えば、以下のような分布として扱います。

$\displaystyle a \sim \mathcal{N}(0,\sigma_a^2), b \sim \mathcal{N}(0,\sigma_b^2)$

これは、 $a$ は平均 $0$ と標準偏差 $\sigma_a$ で、 $b$ は平均 $0$ と標準偏差 $\sigma_b$ で正規分布していることを意味します。

線形回帰モデルには、 $a$ と $b$ 以外にもパラメータがあります。 $y = a + b x + \epsilon , \epsilon \sim \mathcal{N}(0,\sigma^2)$ の $\sigma^2$ です。

この分散 $\sigma^2$ にも何かしら分布を仮定する必要があります。分散なので非負です。そのため、指数分布やガンマ分布、半正規分布など、非負の値を取る分布を仮定します。

分散の場合、指数分布と仮定するケースが多いです。

$\displaystyle \sigma^2 \sim EXP(\lambda)$

ベイズ線形回帰では、このようにパラメータに対し事前に何らかの分布を仮定します。データ量が少ないときには、事前に設定した分布からの影響が大きいですが、学習データが増えれば増えるほどその影響が小さくなります。

データセット

statsmodels経由で取得できるRの以下の、ボストン近郊の住宅価値のデータを利用します。

Housing Values in Suburbs of Boston
https://vincentarelbundock.github.io/Rdatasets/doc/MASS/Boston.html

14変数、506レコードのサンプルデータです。以下の2変数を利用します。1つが目的変数で、もう1つが説明変数です。

目的変数y：ｍedv（持ち家の価格の中央値、千円単位）
説明変数x：lstat（下位ステータスの人の占める割合）

要は、単回帰モデルを構築していきます。

準備

必要なライブラリーを読み込みます。

以下、コードです。

import numpy as np
import pandas as pd
import statsmodels.api as sm

import pymc3 as pm
import arviz as az

from matplotlib import pyplot as plt

# グラフのスタイルとサイズ
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = [12, 9]

必要なデータセットを読み込みます。

以下、コードです。

# データセットの読み込み
dataset = sm.datasets.get_rdataset("Boston", "MASS")

dataset.data #確認

以下、実行結果です。

必要なデータを抜き出し、目的変数yと説明変数xを設定します。

以下、コードです。

y = dataset.data.medv
x = dataset.data.lstat

通常の線形回帰モデル

通常の線形回帰モデルを構築します。

以下、コードです。

# 単回帰モデル
X = sm.add_constant(x) #切片追加

model = sm.OLS(y, X)   #モデルのインスタンス生成
results = model.fit()   #mモデルの学習

results.summary() #確認

以下、実行結果です。

以下の単回帰式が得られました。

$y = 34.5538 - 0.9500 \times x$

グラフ化し見てみます。

以下、コードです。

# 回帰係数と切片の推定値
a = results.params[0] #切片
b = results.params[1] #係数
 
# グラフ化
plt.plot(x, y,"o", label='data') #散布図
plt.plot(x, a + b * x, label='y = a + b * x') #直線
plt.legend()
plt.show()

以下、実行結果です。

ベイズ線形回帰モデル

PyMC3を使いベイズ線形回帰モデルを構築していきます。

　パラメータの事後分布

先ず、パラメータ $a$ 、 $b$ 、 $\sigma^2$ の事後分布を求めてみます。

以下、コードです。

# 各パラメータの事後分布

with pm.Model() as model:
    
    # 事前分布
    a_ = pm.Normal('intercept', a, 10)
    b_ = pm.Normal('slope', b, 10)
    e_ = pm.Exponential('error', 1)
    
    # 予測
    y_ = pm.Normal('observation', a_ + b_ * x, e_, observed=y)
    
    # 事後分布からデータ生成（MCMC）
    trace = pm.sample(return_inferencedata=True)

事後分布のサマリーを見てみます。

以下、コードです。

az.summary(trace)

以下、実行結果です。

グラフで見てみます。

以下、コードです。

az.plot_posterior(trace, point_estimate='mode')

以下、実行結果です。

MCMCのサンプリング状況を見ることもできます。

以下、コードです。

az.plot_trace(trace)

以下、実行結果です。

事後分布から目的変数yのサンプルを生成し、プロットしてみます。

以下、コードです。

# 事後分布から目的変数yのサンプル生成
y_ = pm.sample_posterior_predictive(trace, samples=1000, model=model)

予測値だけ抽出し、y_pred に格納します。

以下、コードです。

y_pred = y_['observation']

y_pred は2次元の配列です。どのような配列なのか見てみます。

以下、コードです。

y_pred.shape

以下、実行結果です。

1,000行×506列の配列です。各 x の予測値 y はどこを見ればいいのでしょうか。

各列が、各 x に該当します。行は、同じx に対する予測 y の値になります。要は、1つの x に対し1,000個の予測値 y があることになります。これが事後分布からのサンプルになります。

各x に対し1,000個の予測値 y の代表値（例：平均値や中央値、最頻値など）を予測値とします。さらに、信頼区間や標準偏差なども計算できます。

ある1つのxに対し、どのような分布になっているのかを見てみます。

ヒストグラムで見てみます。

以下、コードです。

# ヒストグラム表示
plt.hist(y_pred[:,0])

以下、実行結果です。

ついでに、seabornを使ってカーネル密度推定で見てみます。

以下、コードです。

# カーネル密度推定
import seaborn as sns
sns.kdeplot(data=y_pred[:,0])

以下、実行結果です。

1つ1つのxに対し、このように見ていくのは大変ですが、このように分布を出力することができます。

元のデータと、回帰直線、信頼区間をプロットしてみます。

以下、コードです。

# 回帰係数と切片の推定値
a_m = trace.posterior['intercept'].values.mean()
b_m = trace.posterior['slope'].values.mean()

# 信頼区間
x_id = np.argsort(x)
x_ = x[x_id]
ci_95 = az.hdi(y_['observation'], hdi_prob=0.95)[x_id]

# プロット
plt.scatter(x, 
            y, 
            label='data'
           )
plt.plot(x_, 
         a_m + b_m * x_, 
         color = 'b',
         label='y = a + bx'
        )
plt.fill_between(x_, 
                 ci_95[:,0], ci_95[:,1],
                 alpha=0.2, 
                 label='95%HPD'
                )
plt.legend()
plt.show()

以下、実行結果です。

パラメータ $a$ と $b$ の値と予測値yは、分布のサンプルから平均値を求め、その値を利用しています。信頼区間は、予測値yの分布のサンプルから計算しています。

　予測モデル

次に、予測で使えるベイズ線形回帰モデルを構築してみます。外部から与えられた説明変数xに対し予測値を出すモデルです。

以下、コードです。

# 予測モデル
with pm.Model() as predictive_model:
    
    # 事前分布
    a_ = pm.Normal('intercept', a, 10)
    b_ = pm.Normal('slope', b, 10)
    e_ = pm.Exponential('error', 1)
    
    #外部から入力するxのハコ
    x_ = pm.Data('features', x)
    
    # 予測
    y_ = pm.Normal('observation', a_ + b_ * x_, e_, observed=y)
    
    trace = pm.sample()

先ほどと異なり、外から説明変数 $x$ をインプットできるハコのようなものを設定しています（上のコードの中で背景色の異なる箇所）。

新しい説明変数 x の値を作成します。

以下、コードです。

# 予測で利用するx
x_new = np.linspace(0, 40, 400)

0～40の範囲から説明変数用のデータを等間隔で400個作成し、x_new に格納しています。

この x_new をインプットし y を予測してみます。

以下、コードです。

# x_newに対する予測
with predictive_model:
    pm.set_data({'features': x_new})
    posterior = pm.sample_posterior_predictive(trace)

予測値だけ抽出し、y_pred に格納します。

以下、コードです。

# 予測結果（予測値の分布からのサンプル）
y_pred = posterior['observation']

y_pred は2次元の配列です。どのような配列なのか見てみます。

以下、コードです。

# y_predの行×列
y_pred.shape

以下、実行結果です。

4,000行×400列の配列です。

各列が、各 x_new に該当します。行は、同じx_new に対する予測 y の値になります。要は、1つの x_new に対し4,000個の予測値 y があることになります。これが事後分布からのサンプルになります。

各x_new に対し平均値と標準偏差、95%信頼区間の下と上を計算し、データフレームとして整えます。

以下、コードです。

# 各予測の平均値と標準偏差、信頼区間
y_mean = y_pred.mean(axis=0)
y_std = y_pred.std(axis=0)
ci_95 = az.hdi(y_pred, hdi_prob=0.95)

# データフレーム化
pred_result = pd.DataFrame(np.vstack([y_mean,y_std,ci_95[:,0],ci_95[:,1]]).T)
pred_result.columns= ['mean','std','lower','upper']

pred_result #確認

以下、実行結果です。