読者です 読者をやめる 読者になる 読者になる

Obey Your MATHEMATICS.

機械学習関連の純粋数学や実験など

"PyMC3" x "DNN" x "ADVI"で学習させて重みを可視化したら綺麗だったと言う話

こんにちは。


先日今話題沸騰中の ライブラリ Edward でVariational InferenceでBayesian DNNを学習させてみたと言う記事を書きました:

mathetake.hatenablog.com


今回の記事は、Edwardの対抗馬(?)の一つであるPyMC3 を使って

DNNの分類器をADVIでミニバッチ学習させるまでの道のりとその結果を簡単にご報告します。


例のごとく、Higgs粒子のデータセット

mathetake.hatenablog.com

を使います。

PyMCの開発者のブログ記事↓をかなり参考にしました。
Bayesian Deep Learning


コード全文はこの記事の下にGistのリンク張っておきますのでそちらを。

§1 準備

まず諸々のモジュールをimportしてから
前回の記事で作成したデータセットを読み込み、加工します。

import pymc3 as pm
from pymc3 import Model, Normal
import pandas as pd
import theano.tensor as T
import theano
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('white')


#ファイルの読み込み
train_df = pd.read_csv("outlier_removed_high_feature_train.csv", index_col=0, header=0)
test_df = pd.read_csv("high_feature_test.csv", header=0, index_col=0, names=['is_Higgs','m_jj', 'm_jjj', 'm_lv', 'm_jlv', 'm_bb', 'm_wbb', 'm_wwbb'])
train_df = train_df.reindex(np.random.permutation(train_df.index)).reset_index(drop=True)

#トレーニングデータの目的変数をまず-1を0に変換してから抽出
train_df.loc[(train_df.is_Higgs == -1 ),"is_Higgs"] = 0
Y_train = train_df.values[0::,0]

#トレーニングデータの説明変数の正規化
train_df_means = train_df.mean().values
train_df_max = train_df.max().values
train_df_min = train_df.min().values
X_train = train_df.apply(lambda x: (x-x.mean())/(x.max()-x.min()), axis=0).fillna(0).values[0::,1::]

#テストデータの目的変数もまず-1を0に変換してから抽出
test_df.loc[(test_df.is_Higgs == -1 ),"is_Higgs"] = 0
Y_test = test_df.values[0::,0]

#テストデータの説明変数を上の正規化定数を使って正規化
X_test = test_df.values[0::,1::]
for i in range(X_test.shape[1]):
  X_test[0::,i] = (X_test[0::,i] - train_df_means[i+1])/(train_df_max[i+1]-train_df_min[i+1] )

あとは、ネットワークの構成で使う定数を定義しておきます。

total_size = len(Y_train) #トレーニングデータの総数
n_feature = X_train.shape[1] #特徴量の数
n_first_layer_hidden = 20 #第一隠れ層のニューロンの数
n_second_layer_hidden = 15 #第二隠れ層のニューロンの数
mini_batch_size = 10000 #ミニバッチサイズ

§2 ミニバッチの作成

次にミニバッチを返すgeneratorを定義し、ミニバッチを作成します。

# Generator that returns mini-batches in each iteration
def create_minibatch(data, mini_batch_size=50,seed=0):
    rng = np.random.RandomState(seed)   
    while True:
        ixs = rng.randint(len(data), size=mini_batch_size)
        yield data[ixs]
minibatches = zip(
    create_minibatch(X_train, mini_batch_size), 
    create_minibatch(Y_train, mini_batch_size),
)

§3 モデルの構成

準備が出来たのでモデルの構成をします。

ですがモデル構築の前に、
そのInputとtargetになるSharedテンソルを作成しておきます。

ann_input = theano.shared(X_train)
ann_output = theano.shared(Y_train)

ただ学習させるだけであれば、sharedテンソルを用意する必要はありません*1
ですが、今回はテストデータに対する予測値も出したいので、
このX_test,Y_testを放り込むためにsharedテンソルを使います。*2


今回はEdwardの時と同じ4層(7-20-15-1)のネットワークを学習させます。*3
まず事前分布を設計してきます。

with pm.Model() as neural_network:
    #事前分布 on weights (標準正規分布で正則化)
    weights_in_1 = Normal('w_in_1', mu=0, sd=1, shape=(n_feature, n_first_layer_hidden) )
    weights_1_2 = Normal('w_1_2', mu=0, sd=1, shape=(n_first_layer_hidden, n_second_layer_hidden))
    weights_2_out = Normal('w_2_out', mu=0, sd=1, shape=(n_second_layer_hidden,))

    #事前分布 on bias (標準正規分布で正則化)
    bias_on_1 = Normal("b_on_1", mu=0, sd=1,shape=(n_first_layer_hidden,))
    bias_on_2 = Normal("b_on_2", mu=0, sd=1,shape=(n_second_layer_hidden,))

標準正規分布を事前分布として正規化の役目を果たしてもらいます。

次に実際にネットワークを構築します。

with neural_network:
    act_1 = T.tanh(T.dot(ann_input, weights_in_1) + bias_on_1)
    act_2 = T.tanh(T.dot(act_1, weights_1_2) + bias_on_2)
    act_out = T.nnet.sigmoid(T.dot(act_2, weights_2_out))

    # Binary classification -> Bernoulli likelihood
    out = pm.Bernoulli('out', act_out, observed=ann_output)

アウトプットはベルヌーイ分布です。
インプットが上で定義したsharedテンソルann_input
Bernoulliobserved 変数にsharedテンソルann_outputが渡されています。

以上でモデル構築完了です。

§3 ミニバッチADVIのセッティング

学習のセッティングに入ります。

with neural_network:
    minibatch_RVs = [out]
    minibatch_tensors = [ann_input, ann_output]

    #minibatch_ADVIで推定
    v_params = pm.variational.advi_minibatch(n=50000, 
        minibatch_tensors=minibatch_tensors, 
        minibatch_RVs=minibatch_RVs,
        minibatches=minibatches, 
        total_size=total_size, 
        learning_rate=1e-2, epsilon=1.0)

まず、minibatch_RVsで、ミニバッチ学習させる確率変数を格納します。
さらにminibatch_tensors上で定めた確率変数で使われているテンソルの中で、ミニバッチに分割するテンソルを格納します。

この2つと、そのミニバッチ化を実現する先ほど作成したジェネレーターminibatches
pm.variational.advi_minibatch に渡すことでADVIのミニバッチ学習ができるわけです。

§4 学習開始

上で構築したモデルを学習させます。
コードはただ1文で終わりです。
このコードにより、v_paramsで定められた学習が実行され、
学習終了後に事後分布から5000個のサンプルが生成され、traceに格納されます。

with neural_network:
    trace = pm.variational.sample_vp(v_params, draws=5000)

こんな感じのプログレスバーが出てくるはずです↓
f:id:mathetake:20170122155031p:plain

§5 学習結果の確認

学習の過程を確認するには以下のコードで行けます。

plt.plot(v_params.elbo_vals)
plt.ylabel('ELBO')
plt.xlabel('iteration')

f:id:mathetake:20170122155123p:plain

なんだか上手い具合に学習してそうな感じがしますね。


次にパラメータのサマリを確認します。

pm.traceplot(trace);

すると次のように重みの分布が可視化されます。
f:id:mathetake:20170122155630p:plain
f:id:mathetake:20170122155637p:plain
各WeightやBiasの中の各成分が色分けされています。
左側の列は密度推定で可視化された分布図です。右側はサンプリングの回数を横軸に取ったパラメータの分布です。
デルタ関数みたいになっているので、綺麗に最適化(or局所最適化)されていると考えられるでしょう。

綺麗!

ちなみにミニバッチサイズ=2,000で、ADVIのパラメータ更新の回数を1/10にした際のサマリー(一部)は
こんな感じでした↓笑

f:id:mathetake:20170122160738p:plain
全然学習できていません。


さて最後になりましたが、テストデータに対する予測を算出しましょう。

まずsharedテンソルにテストデータを放り込みます。

ann_input.set_value(X_test)
ann_output.set_value(Y_test)

次に事後分布からサンプルを生成します。

#(パラメータの)サンプルを事後分布から生成
ppc = pm.sample_ppc(trace, model=neural_network, samples=50)

最後にppcを使ってAccuracyを算出します。

#アウトプットのサンプルを取り出し、サンプルの平均を取る。
#そして0.5より大きければTrueを取るテストデータサイズのベクトルをpredに格納。
pred = ppc['out'].mean(axis=0) > 0.5

#Accuracyの算出
print('Accuracy = {}%'.format((Y_test == pred).mean() * 100))

得られた結果は
f:id:mathetake:20170122160422p:plain

と言う感じです。
頑張ってハイパーパラメータ調節をして通常の最適化手法で学習した際Accuracyはだいたい70%ほどでしたので、
一発でこの結果が出たと言うのはまずまずと言った感じです。

もっとしっかりと調節すれば良い結果が得られるかもしれません。


楽しい!!!楽しいぞベイズ学習!!



今回の記事は以上です。

PyMCは日本語の詳しい情報がほとんど存在しませんが、
慣れれば簡単にDNNを作れますので、ぜひ一度DNNをベイズ学習させてみてはどうでしょうか?


コード全文↓
PyMC3でDNNをADVIで学習 · GitHub

*1:この場合モデルにデータをそのまま放り込めば良いです。

*2:詳しくは開発者のブログを見て下さい。

*3:全然ディープじゃないやんけ!って方。御尤も。