データサイエンティスト(仮)

元素粒子論博士。今はデータサイエンティスト(仮)。

Pythonでデータ分析:PyStanで線形回帰モデル

導入

ベイズ推定を行うための道具として、マルコフ連鎖モンテカルロMCMC)があります。その派生系であるハミルトニアンモンテカルロ(HMC)をベースにしたソフトウェアとして、Stanというものがよく知られています。

StanはC++ベースのソフトウェアですが、RやPythonを介しても使用することができます。今回は、PythonのインターフェースであるPyStanを使ってベイズ推定を行ってみます。

動機

Stan(特にRStan)は新人時代(2年半くらい前)に触っていたのですが、PJで機械学習の手法を使うことが多くなりあまり使わなくなり…といった感じでしばらく触れていませんでした。最近は、Pythonをベースに手を動かして幅を広げてみようと思っており、その一環としてPyStanに触れるのは良いだろうと思ったのが動機です。というわけで、基本的なことのメモ程度の紹介になると思います。

参考記事

流れを掴んだり、コードを書くに当たって、以下の記事が参考になりました。
qiita.com

インストール

anacondaなどでPythonを導入していると、pipですぐインストールできます。

$ pip install pystan

使用データセット:Boston

久しぶりにボストン近郊の住宅情報のデータを使用します。

import sklearn
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split

Boston = load_boston()
X = Boston.data
y = Boston.target
# 学習用・検証用にデータを分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2)

Stan部分

Stanに渡す部分を以下のように書きます。

# データ(辞書型)
dat = {'N': X_train.shape[0], 'M': X_train.shape[1], 'X': X_train, 'y': y_train}

# Stanコード
model = """
    data {
        int<lower=0> N;
        int<lower=0>M;
        matrix[N, M] X;
        vector[N] y;
    }
    parameters { 
        real beta_0;
        vector[M] beta;
        real<lower=0> sigma;
    }
    model { 
        for (i in 1:N)
            y[i] ~ normal(beta_0 + dot_product(X[i] , beta), sigma);
    }
"""

モデルに与えるデータは辞書型に変更しておきます。

Stanのコードは大まかに以下の3つの部分に分かれています。

  • data
    • データの型などを定義します。配列の長さを明示的に指定(データ数Nや特徴量の数M)するのが地味に大事です。RやPythonを使っていると忘れがちになりますが、StanはC++がベースなので、ちゃんとメモリ領域を人間が指定する必要があります。説明変数が多いときは、matrixで与えるとコードがすっきりします。
  • parameters
    • 推定するパラメータを定義します。ここも説明変数が多いときは、vectorで係数を与えるとコードがすっきりします。
  • model
    • 統計モデルを定義します。今回は目的変数がMEDV(y), 説明変数が他すべてで、残差が正規分布に従うとした場合の正規線形モデルを考えています。
    • dataやparameterでmatrixやvectorで定義した場合は、各レコードに関して係数とデータの内積(dot_product)をとり、for文で回す、とするとすっきりモデル式が書けます。
    •  y = \beta_0 + X^T \beta  + \epsilon, \ \ \ \epsilon \sim N(0, \sigma)

コンパイル

モデルをコンパイルします。結果はstmというオブジェクトにもたせておきます。

%time stm = StanModel(model_code=model)

# 計算時間
CPU times: user 1.11 s, sys: 64 ms, total: 1.18 s
Wall time: 1min 21s

計算時間も一緒に表示させていましが、1分半くらいかかっています(意外とかかるな。。)。

MCMC

サンプリングを実行します。iterationの回数は5000, 初期のiterationは安定的でないことが多いので推定には使わないことが多いのですが、使わないstep(n_warmup)は1000、マルコフ連鎖の数は2(特に大きな理由はないが)としています。

n_itr = 5000
n_warmup = 1000
chains = 2

# サンプリングの実行
%time fit = stm.sampling(data=dat, iter=n_itr, chains=chains, n_jobs=-1, warmup=n_warmup, algorithm="NUTS", verbose=False)

# 計算時間
CPU times: user 556 ms, sys: 142 ms, total: 697 ms
Wall time: 18min 50s

推定に約19分。ちょっと時間かかりますが、ローカルPCだとこんな程度でしょうか。

推定結果

パラメータの推定結果は、fitを出力させると表示できます。

fit

* 出力結果
nference for Stan model: anon_model_0456a22b1cde927056973d34c99a2678.
2 chains, each with iter=5000; warmup=1000; thin=1; 
post-warmup draws per chain=4000, total post-warmup draws=8000.

           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
beta_0    36.89    0.09   5.55  26.09  33.12  36.87  40.63   48.0   4151    1.0
beta[0]   -0.11  3.9e-4   0.03  -0.18  -0.14  -0.11  -0.09  -0.05   8000    1.0
beta[1]    0.05  1.8e-4   0.01   0.02   0.04   0.05   0.06   0.08   6937    1.0
beta[2]    0.07  7.4e-4   0.07  -0.06   0.02   0.07   0.11   0.19   8000    1.0
beta[3]    2.52    0.01    1.0   0.57   1.84   2.52   3.19   4.49   8000    1.0
beta[4]  -17.87    0.06   4.18 -25.99 -20.75 -17.78 -15.09  -9.63   5632    1.0
beta[5]    3.82  6.6e-3   0.47   2.91    3.5   3.83   4.14   4.74   5064    1.0
beta[6] -2.4e-3  1.6e-4   0.01  -0.03  -0.01-2.4e-3 7.7e-3   0.03   8000    1.0
beta[7]   -1.34  2.7e-3   0.22  -1.76  -1.49  -1.35   -1.2  -0.91   6336    1.0
beta[8]    0.32  9.6e-4   0.07   0.17   0.27   0.32   0.37   0.46   5900    1.0
beta[9]   -0.01  5.4e-5 4.2e-3  -0.02  -0.02  -0.01  -0.01-5.7e-3   6138    1.0
beta[10]  -0.98  2.0e-3   0.14  -1.26  -1.08  -0.98  -0.89   -0.7   5230    1.0
beta[11] 8.6e-3  3.2e-5 2.9e-3 2.9e-3 6.7e-3 8.6e-3   0.01   0.01   8000    1.0
beta[12]  -0.51  7.1e-4   0.06  -0.62  -0.55  -0.51  -0.47   -0.4   6223    1.0
sigma      4.61  1.9e-3   0.17    4.3   4.49   4.61   4.72   4.95   8000    1.0
lp__     -817.4    0.05   2.77 -823.7 -819.0 -817.1 -815.4 -812.9   3431    1.0

Samples were drawn using NUTS at Sat Oct 14 13:21:41 2017.
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).

各パラメーターごとに平均、標準誤差、標準偏差、各パーセンタイル等を見ることができます。尤度に関する情報(lp__)もリストの一番下にあります。サンプリング結果の収束を判定する値として、Rhatというものがあり、Rhatが1.1以下だと十分収束していると判断できるようです。今回はすべての係数がRhat1.0となっているので、収束していると判断します。

推定結果の可視化

大雑把にはPyStanのplotメソッドで可視化できます。

fit.plot()

f:id:tekenuko:20171014132852p:plain

説明変数を行列で与えたのが災いして、係数の分布が見づらくなっています。そのため、もう少しわかりやすく可視化します。

# サンプル列を抽出
la  = fit.extract(permuted=True)
# パラメーター名
names = fit.model_pars 
# パラメーターの数
n_param = np.sum([1 if len(x) == 0 else x[0] for x in fit.par_dims])

# プロット
f, axes = plt.subplots(n_param, 2, figsize=(15, 4*n_param))
cnt = 0
for name in names:
    dat = la[name]
    if dat.ndim == 2:
        for j in range(dat.shape[1]):
            d = dat[:,j]
            sns.distplot(d, hist=False, ax=axes[cnt, 0])
            sns.tsplot(d,   alpha=0.8, lw=.5, ax=axes[cnt, 1])
            cnt += 1
    else:
        # Intercept
        sns.distplot(dat, hist=False, ax=axes[cnt, 0])
        sns.tsplot(dat,   alpha=0.8, lw=.5, ax=axes[cnt, 1])
        cnt += 1

name_list = []
for name in names:
    dat = la[name]
    if la[name].ndim == 2:
        for i in range(dat.shape[1]):
            name_list.append("{}{}".format(name,i+1))
    else:
        name_list.append(name)

for i in range(2):
    for j, t in enumerate(name_list):
        axes[j, i].set_title(t)
plt.show()

f:id:tekenuko:20171014134736p:plain:w500

まあ概ね推定できているでしょう。分布的に幾つかは効果がほぼない変数もあるようです。精度や解釈性を向上させたい場合は、そのような変数はモデルから削除することもありますが、ここではこのまま進めます。

精度評価

得られた係数を用いて学習用・検証用データに対して予実比較を行います。今回は、簡単に係数の分布の平均値を推定された係数の値として予測値算出をします。

# 平均値算出
beta0 = la['beta_0'].mean()
beta = la['beta'].mean(axis = 0)
# 予測値算出
y_pred_train = (X_train * beta).sum(axis = 1) + beta0
y_pred_test = (X_test * beta).sum(axis = 1) + beta0
# 学習用、検証用データに関してMSEとR^2を出力
from sklearn.metrics import mean_squared_error, r2_score
print('MSE Train: %.2f, MSE Test: %.2f' % (mean_squared_error(y_train, y_pred_train), mean_squared_error(y_test, y_pred_test)))
print('R^2 Train: %.2f, R^2 Test: %.2f' % (r2_score(y_train, y_pred_train), r2_score(y_test, y_pred_test)))

# 出力結果
MSE Train: 20.41, MSE Test: 28.30
R^2 Train: 0.75, R^2 Test: 0.70

以前自分が書いたsklearnを使った線形回帰モデルの記事では、R^2 Train : 0.763, Test : 0.659となっていました。データの分割の仕方に依存しそうですが、検証用に関しては若干R^2が改善しています。
tekenuko.hatenablog.com

Next Step

今回はボストン近郊の住宅情報のデータセットを用いてベイズ線形回帰モデルを構築しました。次は、データセットとモデル式を変更してより色々なシーンで使えるようにしていこうと思っています。最近KaggleのTitanicのデータをいじりがちなので、特にTitanicでロジスティック回帰なんかが具体的にやろうと思っていることです。