【統計学】stanでロジスティック回帰の実行を割と詳しく解説してみる(w/ Titanic dataset) #Python

StanをつかってTitanicデータをロジスティック回帰してみて、さらに分類の性能評価を少し行ってみるという記事です。

この記事で使う確率的プログラミング言語「Stan」では分布のパラメーターの推定に、ハミルトニアンモンテカルロ法(HMC法)とNUTSという手法が用いられています。厳密には乱数の発生原理が異なるのですが、もう少しシンプルな手法にマルコフ連鎖モンテカルロ法メトロポリス・ヘイスティングス法(MH法)があります。この動作原理について、私@kenmatsu4が書いた

の２点がありますので、よければ参考としてください。やっていることのイメージを付ける意図であればMH法とHMC法は大きく違わないと考えています。（基本的には乱数サンプリングの効率の違い）

環境

OSX Yosemite 10.10.5
Python 2.7
Anaconda 3.18.9
PyStan 2.8.0

コードの全文はGitHubにアップしてあります。

PyStanの導入

色々うまくいかなかったりして悩んだのですが、これがシンプルでうまくいきました。

conda install pystan

ライブラリの準備

import numpy as np
import numpy.random as rd
import pandas as pd

%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import animation as ani
import matplotlib.cm as cm
import seaborn as sns
sns.set(style="whitegrid", palette="muted", color_codes=True)

from tabulate import tabulate
from time import time

import pystan
from pystan import StanModel

from sklearn import cross_validation
from scipy.stats import gaussian_kde

データの準備

【機械学習】モデル評価・指標についてのまとめと実行( w/Titanicデータセット) で使った、数値に変換済みのTitanicデータを使用します。GitHubに用意しましたので、こちらもご利用下さい。

titanic = pd.read_table("https://raw.githubusercontent.com/matsuken92/Qiita_Contents/master/PyStan-Titanic/data/titanic_converted.csv",
              sep=",", header=0)
titanic.head(10)

ID	survived	pclass	sex	age	sibsp	parch	fare	embarked	class	who	adult_male	deck	embark_town	alone
0	0	3	1	22	1	0	7.25	1	3	1	1	0	3	1
1	1	1	0	38	1	0	71.2833	2	1	2	0	3	1	0
2	1	3	0	26	0	0	7.925	1	3	2	0	0	3	0
3	1	1	0	35	1	0	53.1	1	1	2	0	3	3	0
4	0	3	1	35	0	0	8.05	1	3	1	1	0	3	1
5	0	3	1	28	0	0	8.4583	3	3	1	1	0	2	1
6	0	1	1	54	0	0	51.8625	1	1	1	1	5	3	1
7	0	3	1	2	3	1	21.075	1	3	0	0	0	3	0
8	1	3	0	27	0	2	11.1333	1	3	2	0	0	3	0
9	1	2	0	14	1	0	30.0708	2	2	0	0	0	1	0

データを説明変数dataと目的変数targetに分割します。

target = titanic.ix[:, 0]
data = titanic.ix[:, [1,2,3,4,5,6]]

# 訓練データ(80%), テストデータ(20%)に分割する
X_train, X_test, y_train, y_test = cross_validation.train_test_split(data, target, test_size=0.2, random_state=0)
print [d.shape for d in [X_train, X_test, y_train, y_test]]

out

[(712, 6), (179, 6), (712,), (179,)]

つまり、訓練データ$X$は、

のようなデータ数 $N=712$、特徴量数 $M=6$の行列となっています。横ベクトル$x_i$が一人分のデータですね。
6つの特徴量は、
pclass 　　sex　　 age　　 sibsp　　 parch　　 fare
です。

教師データはタイタニックの生存情報で生存=1, 帰らぬ人となった=0 の$N$次元ベクトル

$$
y= (y_1, \cdots, y_N)^{{\rm T}}
$$
です。

ロジスティック回帰について

さてこのようなデータがあったときにロジスティック回帰モデル(wikipedeiaリンク)は下記のように表されます。

\operatorname{logit} (p_{i}) = \ln \left({\frac {p_{i}}{1-p_{i}}}\right)
= \beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i} \\
(i=1, 2, \cdots, N)

$p_i$ は$i$番目のレコードデータの生存確率と見なせます。なぜ生存確率と見なせるか、を数式を少し組み替えて見ていきましょう。まず両辺を$\exp(\cdot)$の中に入れて

$$
\left({\frac {p_{i}}{1-p_{i}}}\right) =
\exp ( \beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i} )
$$

これを$p_i$について解くと

$$
p_i = {1 \over 1 + \exp (\ -(\beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i} )\ ) }
$$

ここで標準シグモイド関数
$$
\sigma(x) = {1 \over 1+ e^{-x}}
$$

を利用します。先ほどの式と形が似ており $x$ を $(\beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i} )$ で置き換えると一致し、

$$
p_i = \sigma(\beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i})
$$

として扱えますね。

グラフにするとこのような形です。

def sigmoid(x):
    return 1 / (1+np.exp(-x))

x_range = np.linspace(-10, 10, 501)

plt.figure(figsize=(9,5))
plt.ylim(-0.1, 1.1)
plt.xlim(-10, 10)

plt.plot([-10,10],[0,0], "k", lw=1)
plt.plot([0,0],[-1,1.5], "k", lw=1)
plt.plot(x_range, sigmoid(x_range), zorder=100)

グラフで見るとわかるようにシグモイド関数$\sigma(x)$は値を(0, 1)の範囲で取るため、この値を確率として解釈して使うことができるのです。

以上より、通常の線形回帰に見える部分、

$$
\beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i}
$$

をベースとしてシグモイド関数をうまく使うことで

\operatorname{logit} (p_{i}) = \ln \left({\frac {p_{i}}{1-p_{i}}}\right)
= \beta + \beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i} \\
(i=1, 2, \cdots, N)

の形に持って行き、値を確率と解釈して使うのがロジスティック回帰と言えます。

この式のパラメーター$\beta, \beta_1, \cdots, \beta_M$ を今回Stanを使って推定してみようということで、Stanが登場するわけです。

Stanでの実行

Stanに渡すために、訓練データを詰めたDictionaryオブジェクトを作成します。

dat = {'N': X_train.shape[0], 'M': X_train.shape[1], 'X': X_train, 'y': y_train}

code = """
        data {
            int<lower=0> N;
            int<lower=0> M;
            matrix[N, M] X;
            int<lower=0, upper=1> y[N];
        } parameters {
            real beta0;
            vector[M] beta; 
        } model {
            for (i in 1:N)
                y[i] ~ bernoulli(inv_logit (beta0 + dot_product(X[i] , beta)));
        }
        """

※ stanコードはこちらの記事を参考にさせていただきました。

Stanのコードは幾つかのブロックに分かれていますが今回は

data
parameters
model

の3つのパートに分かれています。

1つずつ見ていきます。

dataブロック

data {
    int<lower=0> N;
    int<lower=0> M;
    matrix[N, M] X;
    int<lower=0, upper=1> y[N];
}

ここは、Python側から値を詰め込む先の、いわゆるプレースホルダーのような役割です。Python側で

dat = {'N': X_train.shape[0], 'M': X_train.shape[1], 'X': X_train, 'y': y_train}

のように値を代入していたので、それがStan側で使われます。
今回、

N: データ数
M: 特徴量数
X: NxM行列の訓練データ
y: N次元の教師データ。上限値1, 下限値0という制限を入れている。

のように各変数を使っています。

parametersブロック

parameters {
    real beta0;
    vector[M] beta; 
}

parametersブロックでは、確率分布のパラメーターなど、推定したいターゲットの変数を宣言するブロックです。今回は回帰パラメーター（切片 $\beta$と回帰係数 $\beta_i \ (i=1,2,\cdots,M)$)を推定したいのでここで宣言します。

modelブロック

ここがキモですね。統計モデルの宣言です。

model {
    for (i in 1:N)
        y[i] ~ bernoulli(inv_logit (beta0 + dot_product(X[i] , beta)));
}

Stanの関数が幾つかあるのでまずそれを紹介します。

real dot_product(vector x, vector y)
Stan Manual v2.9 p358
　内積を計算します。ここでは数式の $\beta_{1}x_{1,i} + \cdots + \beta_{M}x_{M,i}$ の部分に該当し、各レコードとパラメータをかけて足し合わせています。

real inv_logit(real y)
Stan Manual v2.9 p338
下記の数式で表される計算をする関数です。つまり標準シグモイド関数ですね。
$$
\operatorname{inv \_ logit}(y) = {1 \over 1 + \exp(−y)}
$$

y ~ bernoulli(theta)
Stan Manual v2.9 p384
~で表されるこの式はSampling Statementsと呼ばれ、実際は下記の式、
increment_log_prob(bernoulli_log(y, theta));
が呼び出され実行されています。ベイズ推定による事後確率を計算するので、観測済みデータとしてyを、推定したいパラメーターとしてthetaを引数に入れています。thetaは先ほど説明をした説明変数とパラメータ$\beta_i$を線形結合して標準シグモイド関数の引数に設定していたinv_logit (beta0 + dot_product(X[i] , beta))の部分ですね。
まとめると、この１行は対数尤度を計算しているパートであり、bernoulli分布の確率関数が

$$
B(y_i|\theta) = \theta^{y_i}(1-\theta)^{1-y_i}
$$

であり、その対数をとったものは、

$$
\ln(B(y_i|\theta)) = y_i \ln (\theta) + (1-y_i) \ln(1-\theta)
$$

さらに、データ数 $N$の同時確率なので、対数をとると足し算となり、

$$
\sum_{i=1}^N\ln(B(y_i|\theta)) = \sum_{i=1}^N\ \{ y_i \ln (\theta) + (1-y_i) \ln(1-\theta) \}
$$

です。この足し算部分がincrement_log_probで処理されている部分で、for文で回しながら積み上げで足されています(なのでincrement)。

また、パラメーターbeta0, betaに対して事前分布が明示的に設定されていませんが、この場合Stanによって無情報事前分布（範囲(−∞, ∞)の一様分布）が設定されています。

以上のStanコードから事後分布に比例する値(尤度×事前分布)の用意ができました。これを使ってHMC,NUTSによる乱数生成を行います。

参考：Stan Manual v2.9
　　　26.3. Sampling Statements
　　　6.2. Priors for Coefficients and Scales

モデルのコンパイル

StanはそのコードをC++のコードに変換してコンパイル、それを実行することで高速化を行っています。毎回コードのコンパイルをするのは非効率なので、コードが一度決まればコンパイルしてPythonのオブジェクトとして持っていた方が何かと便利です。

# Build Stan model
%time stm = StanModel(model_code=code)

out

CPU times: user 1.36 s, sys: 113 ms, total: 1.48 s
Wall time: 15.4 s

モデルのビルドが完了したら、サンプリングの実行です。ここが一番時間のかかるところですね。
サンプリング数に関する設定は下記の通りです。

サンプリング数: 3000
ウォームアップ数¹(バーンイン数): 200
chain数: 200

(3000-200)*2 = 5600個のサンプル数となるまでサンプリングを行います。

実行の設定としては
* 並列実行数: -1を設定することでマシンの可能最大並列数で実行します。
* サンプリングのアルゴリズムにはNUTSを使用します。

n_itr = 3000
n_warmup = 200
chains = 2

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

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

サンプリングの実行には５分ほどかかりました。

out

CPU times: user 156 ms, sys: 59.1 ms, total: 215 ms
Wall time: 5min 24s

サンプリング結果のサマリーを表示します。

print fit

各パラメーターごとに平均、標準誤差、標準偏差、各パーセンタイル等を見ることができます。尤度に関する情報も最下部に示されています。サンプリング結果の収束を判定する値Rhatが1.0なので、十分分布が収束していることがわかります。(一般にRhatが1.1以下であれば十分収束していると判断できるとのことです² )

out

Inference for Stan model: anon_model_b12e3f2368679a0c562b9f74618b2f82.
2 chains, each with iter=3000; warmup=200; thin=1; 
post-warmup draws per chain=2800, total post-warmup draws=5600.

          mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
beta0     5.02  8.0e-3    0.6   3.89   4.61   5.01   5.42   6.23 5600.0    1.0
beta[0]  -1.04  2.1e-3   0.15  -1.35  -1.15  -1.04  -0.93  -0.75 5600.0    1.0
beta[1]  -2.77  3.0e-3   0.23  -3.22  -2.92  -2.76  -2.61  -2.33 5600.0    1.0
beta[2]  -0.04  1.2e-4 8.9e-3  -0.06  -0.05  -0.04  -0.04  -0.03 5600.0    1.0
beta[3]  -0.41  1.6e-3   0.12  -0.66  -0.49  -0.41  -0.33  -0.18 5600.0    1.0
beta[4]  -0.08  1.8e-3   0.14  -0.35  -0.17  -0.08   0.02   0.18 5600.0    1.0
beta[5] 2.5e-3  3.4e-5 2.5e-3-2.3e-3 7.0e-4 2.4e-3 4.1e-3 7.7e-3 5600.0    1.0
lp__    -322.4    0.02   1.86 -326.8 -323.4 -322.0 -321.0 -319.7 5600.0    1.0

Samples were drawn using NUTS(diag_e) at Sun Dec 13 02:53:45 2015.
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).

各パラメーターの分布の期待値をEAP(Expected A Posteriori : 事後期待値)推定量と言います。そのEAP推定量は単にサンプル結果の平均をとることで得られます。簡単ですね

# 各パラメーターのEAP推定量リスト抽出
mean_list = np.array(fit.summary()['summary'])[:,0]

次に、MAP(Maximum A Posteriori : 最大事後確率)推定量を求めます。これはちょっと工夫しまして、カーネル密度推定を行ってその最大値をとる値を持ってきて算出しました。（もっと簡単な方法を知っている方がいれば教えて下さいm(_ _)m ）

# 各パラメーターのMAP推定量のリスト生成
ddd = la['beta0']
def find_MAP(data):
    kde = gaussian_kde(data)
    x_range = np.linspace(min(data), max(data), 501)
    eval_kde = kde.evaluate(x_range)
    #plt.plot(x_range, eval_kde)
    return x_range[np.argmax(eval_kde)]

MAP_list = []
MAP_list.append(find_MAP(ddd))
for i in range(n_param-1):
    MAP_list.append(find_MAP(la['beta'][:,i]))

最後に各パラメーターのサンプリング結果をグラフで可視化します。
左側はサンプリング結果をカーネル密度にしたもの、右側はサンプリングした順番にその乱数の値をプロットしたtrace plotと呼ばれるものです。

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, rug=True, 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, rug=True, ax=axes[cnt, 0])
        sns.tsplot(dat,   alpha=0.8, lw=.5, ax=axes[cnt, 1])
        cnt += 1

name_list = []
for name in names:
    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()

同じく、尤度に関してグラフを書いたものです。

# Likelihood
f, axes = plt.subplots(1, 2, figsize=(15, 4))

sns.distplot(la['lp__'], hist=False, rug=True, ax=axes[0])
sns.tsplot(la['lp__'],   alpha=0.8, lw=.5, ax=axes[1])
axes[0].set_title("Histgram of likelihood.")
axes[1].set_title("Traceplot of likelihood.")
plt.show()

最後に推定結果の性能を見てみます。最初にデータをホールドアウト法で8:2に分けていたので、再代入した時の性能と、汎化性能をみてみます。使用するパラメーターとして、EAP推定量、MAP推定量ともに使って比べてみます。

# 推定したパラメータ値を代入して確率p_iを算出するための関数。
def logis(x, beta):
    assert len(beta) == 7
    assert len(beta) == 7
    if type(x) != 'np.array':
        x = np.array(x)
    tmp = [1]
    tmp.extend(x)
    x = tmp
    return (1+np.exp(-np.dot(x, beta)))**(-1)

# 設定された閾値で0 or 1を判定。デフォルト閾値は0.5
def check_accuracy(data, target, param, threshold = 0.5):
    ans_list = []
    for i in range(len(data)):
        idx = data.index[i]
        res = logis(data.ix[idx], param)
        ans = 1 if res > threshold else 0
        ans_list.append(ans == target.ix[idx])

    return np.mean(ans_list)


param = mean_list[0:7]

# 再代入した時の正答率
print u"[train][EAP] Accuracy:{0:.5f}".format(check_accuracy(X_train, y_train, param))
print u"[train][MAP] Accuracy:{0:.5f}".format(check_accuracy(X_train, y_train, MAP_list))

# テストデータの正答率で汎化性能を見る
print "[test][EAP] Accuracy:{0:.5f}".format(check_accuracy(X_test, y_test, param))
print "[test][MAP] Accuracy:{0:.5f}".format(check_accuracy(X_test, y_test, MAP_list))

out

[train][EAP] Accuracy:0.79073
[train][MAP] Accuracy:0.80618
[test][EAP] Accuracy:0.81006
[test][MAP] Accuracy:0.79888

再代入時はMAP推定量の方が正答率が高いですが、汎化性能的にはEAP推定量の法がよさそうに見えます。これがたまたまなのか、そのような傾向があるのか興味があるところです。

参考

http://aaaazzzz036.hatenablog.com/entry/2013/11/06/225417
　 → Stanコードを参考とさせていただきました。

Stanマニュアル (Modeling Language User’s Guide and Reference Manual, v2.9.0)
https://github.com/stan-dev/stan/releases/download/v2.9.0/stan-reference-2.9.0.pdf

BDA3(Bayesian Data Analysis 3rd edition p282の注釈)によると、バーンインというのは誤解を生む場合があるので、ウォームアップと呼ぶことを推奨しているようです。 ↩
参考：MCMCの収束診断におけるRhatの具体的な値についての引用メモ ↩