Memo：Gluonの解説やコード紹介（海外）
導入
2017年10月12日（現地時間）に、MicrosoftとAWSがGluonというDeep Learningのライブラリを公開しました。
日本語の解説記事があまり見当たらなかったので、簡単なところは自分で試してみるなどし、いくつか記事にもしました。色々調べていたところ、Documentや国際会議などではすでにGluonじたいやコードの紹介があることがわかってきました。こういった情報を整理しておくのは後々の自分にとって有用だろうと思ったので、メモがてら残しておきます。
GithubのLecture Note
ICML2017, KDD2017での発表者のGithubにMXNetも含めたLecture Noteがあります。形式はJupyter Notebookです。MXNetのメソッドの使い方から、深層強化学習まで、広い範囲を扱っています。これを最初に見つけてたらどんなに楽だったか。
github.com
また、上記のLectureに対応するWeb Pageもあります。
GluonでDeep Learning：CNNを組んでみる
導入
前回、MicrosoftとAWSが公開したライブラリであるGluonの紹介をしました。
tekenuko.hatenablog.com
前回紹介したのは、Tutorialの多層パーセプトロン（MLP）でしたが、Gluonは他のネットワークもサポートしています。今回は、畳み込みニューラルネットワーク（Convolutional Neural Network, CNN）の例を紹介しようと思います。
参考
参考記事などはまだ見当たらなかったので、四苦八苦してたのですが、MXNetのDocumentを探してみると、使用例がありました。
Convolutional Neural Networks in gluon — The Straight Dope 0.1 documentation
概ねこの使用例を踏襲して記述していけばよさそうです。
データセット：MNIST
前回と同じく、MNISTを使います。いったん、バッチサイズや出力層の次元数を変数と読み込んだデータ（MXNetのNDArray形式のデータ）を変換する関数を定義しておきます。nd.transpose()でうまく転置するのがポイントで、自分はこの処理をやってなくてエラーが出て対処に困りました（汗）。
import mxnet as mx from mxnet import nd, autograd from mxnet import gluon import numpy as np batch_size = 64 num_outputs = 10 def transform(data, label): return nd.transpose(data.astype(np.float32), (2,0,1))/255, label.astype(np.float32)
MNISTデータはMXNetのサイトからロードします。
train_data = mx.gluon.data.DataLoader(mx.gluon.data.vision.MNIST(train=True, transform=transform), batch_size, shuffle=True) test_data = mx.gluon.data.DataLoader(mx.gluon.data.vision.MNIST(train=False, transform=transform), batch_size, shuffle=False)
ネットワーク構築
Gluonでは、Sequential()というメソッドを定義し、そこに積み木のように層を足していく感覚でネットワークのモデルを構築できます。今回は畳み込み層とプーリング層を2回挟み、その後全結合層を入れて最後に出力層、といった構成にしました。畳み込み層→プーリング層の組み合わせは、ある程度微小変換に対してロバストな特徴をデータから抽出するのによく用いられます。その後、全結合層を入れることで、ざっくりそういった抽出した特徴をどういった重みで組み合わせるのがよいかを表現しています。
# 全結合層の次元 num_fc = 512 net = gluon.nn.Sequential() with net.name_scope(): net.add(gluon.nn.Conv2D(channels=20, kernel_size=5, activation = 'relu')) net.add(gluon.nn.MaxPool2D(pool_size=2, strides=2)) net.add(gluon.nn.Conv2D(channels=50, kernel_size=5, activation = 'relu')) net.add(gluon.nn.MaxPool2D(pool_size=2, strides=2)) net.add(gluon.nn.Dense(num_fc, activation = 'relu')) # 出力層の次元は10 net.add(gluon.nn.Dense(num_outputs))
畳み込み層、プーリング層を今回は入れています。どんな引数が使えるかは、以下のページを見るとわかります。
github.com
- Conv2D
- 2次元の畳み込みを行う層です。フィルタの数（channel）やサイズ（kernel）が最低限指定する引数になります。
- 他、画像データに対するフィルタのスライドのさせ方（strides）や、入出力の画像のサイズを変えないための処方（padding）のやり方などを指定できます。
- ここはDocumentを見て、デフォルト以外を試したかったら変更するとよいです。
- MaxPool2D
- 最大プーリングを行う層です。特定のWindowの値を最大値で代用し、集約します。
- 何も指定しないとデフォルト（2×2のWindowで最大プーリング、ストライドはNoneなど）
- pool_sizeは数値 or list/tupleで指定できます。今回はあえて数値で指定しています。
学習の設定
前回とほぼ同じ設定です。まず、CPUを使用することを宣言します。
ctx = mx.cpu()
初期値はXavierの初期値というものを今回は用いています。目的関数や最適化手法は変更していません。
net.collect_params().initialize(mx.init.Xavier(magnitude=2.24), ctx=ctx) softmax_cross_entropy = gluon.loss.SoftmaxCrossEntropyLoss() trainer = gluon.Trainer(net.collect_params(), 'sgd', {'learning_rate': .1})
初期化はDocumentがやっていたことを変えていないだけです（笑）。Xavierの初期化は、乱数で重みに初期値を与える際に、層のノードの数で重み付けを行う方法になっています。重みのばらつきをノードの数を考慮して均一化し、学習をうまく進めるよう工夫しているようです。解説に関しては、以下のようなページがあります。
qiita.com
学習
大まかな設定は前回と同じです。今回は、学習データとテストデータのAccuracyを出力できるようにしました。データとモデルを指定するとAccuracyを出力する関数を予め作っておきます。
def evaluate_accuracy(data_iterator, net): acc = mx.metric.Accuracy() for i, (data, label) in enumerate(data_iterator): data = data.as_in_context(ctx) label = label.as_in_context(ctx) output = net(data) predictions = nd.argmax(output, axis=1) acc.update(preds=predictions, labels=label) return acc.get()[1]
epoch数は10で学習を進めます。epochごとにlossと学習データ、テストデータのAccuracyを出力します（lossは移動平均も考えていたりしてますが、Documentに書いてあることをそのまま書いてあるだけです（笑））。
epochs = 10 smoothing_constant = .01 for e in range(epochs): for i, (data, label) in enumerate(train_data): data = data.as_in_context(ctx) label = label.as_in_context(ctx) with autograd.record(): output = net(data) loss = softmax_cross_entropy(output, label) loss.backward() trainer.step(data.shape[0]) ########################## # 損失関数の移動平均 ########################## curr_loss = nd.mean(loss).asscalar() moving_loss = (curr_loss if ((i == 0) and (e == 0)) else (1 - smoothing_constant) * moving_loss + (smoothing_constant) * curr_loss) test_accuracy = evaluate_accuracy(test_data, net) train_accuracy = evaluate_accuracy(train_data, net) print("Epoch %s. Loss: %s, Train_acc %s, Test_acc %s" % (e, moving_loss, train_accuracy, test_accuracy)) # 出力 Epoch 0. Loss: 0.0838328978149, Train_acc 0.980033333333, Test_acc 0.9826 Epoch 1. Loss: 0.0540922844512, Train_acc 0.985, Test_acc 0.9856 Epoch 2. Loss: 0.0414136623197, Train_acc 0.991266666667, Test_acc 0.9897 Epoch 3. Loss: 0.0309304529114, Train_acc 0.989783333333, Test_acc 0.9868 Epoch 4. Loss: 0.0209868983092, Train_acc 0.995583333333, Test_acc 0.9917 Epoch 5. Loss: 0.0190144998394, Train_acc 0.988766666667, Test_acc 0.9848 Epoch 6. Loss: 0.0142381958556, Train_acc 0.997033333333, Test_acc 0.9915 Epoch 7. Loss: 0.0122593285998, Train_acc 0.997483333333, Test_acc 0.9912 Epoch 8. Loss: 0.0127596473961, Train_acc 0.997866666667, Test_acc 0.9915 Epoch 9. Loss: 0.00813980522647, Train_acc 0.9982, Test_acc 0.9913
ローカルPCだと計算に20分くらいかかります。ですが、学習、テストデータともにAccuracyが99%を超えています。感覚的には、中間層1層のMLPだとAccuracyが90%前半程度になので、CNNにしたことで大幅な精度向上になっていると期待できます。
別のネットワーク例
自作でCNNを作っていたとき、どうしてもエラーが解消できなくてネットサーフィンしてたときに見つけた例を紹介しておきます。
num_fc = 512 net = gluon.nn.Sequential() with net.name_scope(): net.add(gluon.nn.Conv2D(channels=20, kernel_size=5)) net.add(gluon.nn.BatchNorm(axis=1, center=True, scale=True)) net.add(gluon.nn.Activation(activation='relu')) net.add(gluon.nn.MaxPool2D(pool_size=2, strides=2)) net.add(gluon.nn.Conv2D(channels=50, kernel_size=5)) net.add(gluon.nn.BatchNorm(axis=1, center=True, scale=True)) net.add(gluon.nn.Activation(activation='relu')) net.add(gluon.nn.MaxPool2D(pool_size=2, strides=2)) net.add(gluon.nn.Flatten()) net.add(gluon.nn.Dense(num_fc)) net.add(gluon.nn.BatchNorm(axis=1, center=True, scale=True)) net.add(gluon.nn.Activation(activation='relu')) net.add(gluon.nn.Dense(num_outputs))
途中でBatchNormalization（ざっくり学習を効率的にすすめるための工夫です）を施したりしています。上で自分が作ったネットワークは、スタンダードなCNNの要素に絞りたかったので、細かい工夫はカットしてネットワークを構築しました。
このネットワークを学習させると、Accuracyは以下のようになります。
# 学習と検証結果 Epoch 0. Loss: 0.0459737773102, Train_acc 0.992066666667, Test_acc 0.9892 Epoch 1. Loss: 0.0292876055092, Train_acc 0.992816666667, Test_acc 0.9882 Epoch 2. Loss: 0.0204365993603, Train_acc 0.995516666667, Test_acc 0.9908 Epoch 3. Loss: 0.015462121763, Train_acc 0.998183333333, Test_acc 0.9926 Epoch 4. Loss: 0.0106946259829, Train_acc 0.999033333333, Test_acc 0.9926 Epoch 5. Loss: 0.00973990939113, Train_acc 0.999466666667, Test_acc 0.9926 Epoch 6. Loss: 0.00662801339947, Train_acc 0.999466666667, Test_acc 0.993 Epoch 7. Loss: 0.00540244358498, Train_acc 0.9998, Test_acc 0.9927 Epoch 8. Loss: 0.00381258537248, Train_acc 0.999916666667, Test_acc 0.993 Epoch 9. Loss: 0.00272173117527, Train_acc 0.99995, Test_acc 0.9934
工夫しているだけあって、普通にCNN組むよりも性能が上がるようですね。
Next Step
今回は、Gluonで畳み込みニューラルネットワーク（Convolutional Neural Network, CNN）の例を紹介しました。ネットワークの構築は、CNNの場合も同様に非常に簡単に行うことができます。次回以降は、異なるデータを使う、CNN以外のネットワークを試す、などを検討しようと思います。
GluonでDeep Learning：Tutorialを眺めてみる
導入
2017年10月12日（現地時間）に、MicrosoftとAWSがGluonというDeep Learningのライブラリを公開しました。
www.itmedia.co.jp
Gluonとは、自然界の基本的な相互作用の一つ「強い相互作用」を伝える素粒子のことです。glue（のり）にちなんでのりのように強く粒子を結びつけるといった感覚で名付けられました。ちなみに、ゲージ群でいうとSU(3)の随伴表現にアサインされます（）。
立ち位置としては、TheanoやTensorflowに対するKerasに近い印象です。AWSがサポートしてるDeep LearningのフレームワークであるMXNetをサポートしており、コーディングを単純化・簡単化するようなものとなっています。Kerasと異なる点は、MXNetの一部のようになっている点でしょうか。Gluonは、インポートするのはあくまでMXNetで、その中のライブラリとしてGluonがある、といった感じになっています（触ってみた印象です）。
現在、Github上でDocumentやTutorialが公開されています。
github.com
こちらを見れば雰囲気はわかると思いますが、日本での紹介がまだ見当たらなかったので、眺めるついでに紹介をしたいと思います。
インストール方法
anacondaなどでPythonを導入していると、pipですぐインストールできます。ただし、インストールするのはmxnetです。
$ pip install mxnet
おそらく問題なくインストールできると思います。
Tutorial：MNIST
Githubのページには、MNISTで多層パーセプトロン（MLP）の実装例が載っています。
github.com
以下、Tutorialを眺めていきます。
ライブラリのインポート
mxnet、および必要なライブラリをインポートします。Gluonに対応するライブラリはgluonです。numpyは、Deep Learningのフレームワークに入力するデータの型をnp.XX32にする、とかいったデータの扱いまわりで必要なので、一緒にインポートしています。
import mxnet as mx from mxnet import gluon, autograd, ndarray import numpy as np
使用データセット：MNIST
手書き数字の画像と数字のラベルがセットになっているオープンなデータセットであるMNISTを使います。手書き数字の画像は28×28ピクセルで、各ピクセルに0から255の値が入っている2次元配列です（要はグレースケール）。各ピクセルの値を最大1になるように正規化の処理をし、型もnp.float32とします。このような処理をしないと、フレームワークによっては学習がうまく進まなかったりエラーになったりします。また、後で処理しますが、今回はMLPを想定しているので、画像データは28×28=784次元のベクトルデータとします。
Tutorialでは、MNISTのデータをmxnetのサイトからダウンロードしてきています。
train_data = mx.gluon.data.DataLoader(mx.gluon.data.vision.MNIST(train=True, transform=lambda data, label: (data.astype(np.float32)/255, label)), batch_size=32, shuffle=True) test_data = mx.gluon.data.DataLoader(mx.gluon.data.vision.MNIST(train=False, transform=lambda data, label: (data.astype(np.float32)/255, label)), batch_size=32, shuffle=False)
mx.gluon.data.DataLoaderですが、ロードしてきた（or 読み込んだ）データに対して、ミニバッチに分割してくれるメソッドです。ミニバッチのサイズや、シャッフルしてミニバッチにデータをアサインするかなどを指定できます。
trainデータは60000万枚、testデータは10000枚の画像データになっています。
ネットワーク構築
この部分がGluonの売りの一つになると考えられます。Gluonでは、Sequential()というメソッドを定義し、そこに積み木のように層を足していく感覚でネットワークのモデルを構築できます。KerasやCognitive Toolkitと似ていますね。
# モデルの初期化：Sequentialメソッドをインスタンス化 net = gluon.nn.Sequential() # addメソッドで層を追加する with net.name_scope(): net.add(gluon.nn.Dense(128, activation="relu")) net.add(gluon.nn.Dense(64, activation="relu")) net.add(gluon.nn.Dense(10)) # 出力層
Decse()は全結合層を足すメソッドで、ここでは引数としてユニットの数と活性化関数を指定しています。1層目と2層目のユニットの数はそれぞれ128と64で、活性化関数は両方ともReLUです。出力層に関しては、0~9の10種類の数字のクラス分類のタスクを考えるので、出力層の数は10とします。
出力層に関してコメントですが、活性化関数は学習の際に損失関数とセットで指定する思想のようです。ここはCognitive Toolkitと似ています。この思想のキモチは、モデルは「学習されるパラメータがある部分」ということなんだろうと個人的には思っています。クラス分類をする際のソフトマックス層は学習するパラメータはないので、ネットワークの出力を目的関数へ変換する際の一部分という位置づけになっているのだと思います。
学習の設定
ここはかなりCognitive Toolkit的です。パラメータの初期化、損失関数、最適化手法などを決め、それをTrainerというメソッド（どういう学習をするかのスケジューリングを司るもの）に放り込みます。
# パラメータ初期化：標準偏差0.05の正規乱数を利用 net.collect_params().initialize(mx.init.Normal(sigma=0.05)) # 損失関数はクロスエントロピー softmax_cross_entropy = gluon.loss.SoftmaxCrossEntropyLoss() # 最適化手法はSGDで、学習率は0.1 trainer = gluon.Trainer(net.collect_params(), 'sgd', {'learning_rate': .1})
これでネットワークのパラメータを学習する準備ができました。
学習
先程までの設定に従って学習していきます。何回データ学習に使うかの回数であるepoch数は10にしています。前にコメントしましたが、入力はミニバッチごとに投入、入力の形式はreshapeにより1次元のベクトルに変換しています。誤差逆伝搬に関係する微分操作などは、autogradが司っています。
# エポック数は10 epochs = 10 for e in range(epochs): for i, (data, label) in enumerate(train_data): data = data.as_in_context(mx.cpu()).reshape((-1, 784)) label = label.as_in_context(mx.cpu()) with autograd.record(): # 誤差逆伝搬をやる output = net(data) loss = softmax_cross_entropy(output, label) loss.backward() trainer.step(data.shape[0]) curr_loss = ndarray.mean(loss).asscalar() print("Epoch {}. Current Loss: {}.".format(e, curr_loss)) # 出力 Epoch 0. Current Loss: 0.08900555968284607. Epoch 1. Current Loss: 0.2502848505973816. Epoch 2. Current Loss: 0.2667727470397949. Epoch 3. Current Loss: 0.017404090613126755. Epoch 4. Current Loss: 0.006391774397343397. Epoch 5. Current Loss: 0.07984966039657593. Epoch 6. Current Loss: 0.009250413626432419. Epoch 7. Current Loss: 0.062348250299692154. Epoch 8. Current Loss: 0.06544484943151474. Epoch 9. Current Loss: 0.0023109368048608303.
実行すると、上のような感じでもりもり損失関数が小さくなっていきます。
その他、Gluonでできそうなこと
以上がTutorialになります。Tutorialだけだと予測も予実比較もしていないので、ネットワークの性能が知りたい場合は適宜追加で書く必要があります。
Documentを見る限りだと、MLP以外にもネットワークのモデルはサポートしているようです（今の時代当たり前だと思いますが）。ざっと見てみると
- 畳み込みニューラルネットワーク（CNN）
- 1, 2, 3次元の畳み込み層
- プーリング（Max, Avg, GlobalMax）
- 他テクニック
- BatchNormalizationやDropout
- ActivationにLeakyReLUがある
- 再帰型ニューラルネットワーク（RNN）
- RNN
- LSTM
- GRU
- Cell
- RNNCell, GRUCell, RecurrentCell, SequentialRNNCell, BidirectionalCell, DropoutCell, ZoneoutCell, ResidualCell
など、代表的なネットワークは作れる状態のようです。
Tutorialの内容の通り、比較的簡単にDeep Learningを試せるようなライブラリとなっています。また、今はMXNetだけですが、今後はCognitive ToolkitやCaffe2もサポートする予定とのことです。というわけで、複数のフレームワークを横断して使いやすいものになることが期待されるため、これからDeep Learningをやりたいと考えている人は検討するとよいのかもしれません。
Pythonでデータ分析：PyStanで線形回帰モデル
導入
ベイズ推定を行うための道具として、マルコフ連鎖モンテカルロ（MCMC）があります。その派生系であるハミルトニアンモンテカルロ（HMC）をベースにしたソフトウェアとして、Stanというものがよく知られています。
StanはC++ベースのソフトウェアですが、RやPythonを介しても使用することができます。今回は、PythonのインターフェースであるPyStanを使ってベイズ推定を行ってみます。
動機
Stan（特にRStan）は新人時代（2年半くらい前）に触っていたのですが、PJで機械学習の手法を使うことが多くなりあまり使わなくなり…といった感じでしばらく触れていませんでした。最近は、Pythonをベースに手を動かして幅を広げてみようと思っており、その一環としてPyStanに触れるのは良いだろうと思ったのが動機です。というわけで、基本的なことのメモ程度の紹介になると思います。
参考記事
流れを掴んだり、コードを書くに当たって、以下の記事が参考になりました。
qiita.com
使用データセット：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
- parameters
- 推定するパラメータを定義します。ここも説明変数が多いときは、vectorで係数を与えるとコードがすっきりします。
- model
コンパイル
モデルをコンパイルします。結果は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()
説明変数を行列で与えたのが災いして、係数の分布が見づらくなっています。そのため、もう少しわかりやすく可視化します。
# サンプル列を抽出 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()
まあ概ね推定できているでしょう。分布的に幾つかは効果がほぼない変数もあるようです。精度や解釈性を向上させたい場合は、そのような変数はモデルから削除することもありますが、ここではこのまま進めます。
精度評価
得られた係数を用いて学習用・検証用データに対して予実比較を行います。今回は、簡単に係数の分布の平均値を推定された係数の値として予測値算出をします。
# 平均値算出 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
Pythonでデータ分析：Catboost
導入
2017年7月に、ロシアのGoogleと言われている（らしい）Yandex社から、Catboostと呼ばれるGradient Boostingの機械学習ライブラリが公開されています。
ここ何ヶ月か調整さんになっていて分析から遠ざかりがちになりやすくなっていたのですが、手を動かしたい欲求が我慢できなくなってきたので、いい機会だと思って触ってみることにしました。
使用データセット：Titanic
Titanic号の乗客の情報および生存・死亡のデータを用いて、どの程度の精度のモデルができるかを試してみます。
データセットは、Kaggleのページからダウンロードしました。
Titanic: Machine Learning from Disaster | Kaggle
train.csvとtest.csvの2つを利用します。train.csvはSurvived（生存なら1、死亡なら0）のカラムがあるデータ、test.csvはSurvivedのカラムが無いデータになっています。train.csvを用いてモデルを作成し、test.csvからSurvivedの予測値を出し、Kaggleにsubmitすることで、モデルの精度(今回はAccuracy)を見ます。
train.csvのデータはdf_raw、test.csvのデータはdf_raw_testという名前で読み込んでおきます。
データ加工（笑）
Titanicのデータは、いくつか欠損があるカラムがあるのですが、今回は雑に-999で埋めます。また、乗客の名前などの非構造データも幾つかあるのですが、それらは予めインデックスを保持しておき、モデリングの際によしなにカテゴリカル変数として処理してもらうことにします（Catboostはそういった機能があったので、挙動確認も込めてあえて特徴量加工を頑張らないことにしました）。
# pandasやnumpyはインポート済み # 欠損値補完 df_raw.fillna(-999, inplace=True) # 説明変数・目的変数 X = df_raw.drop('Survived', axis = 1) y = df_raw['Survived'] # カテゴリカル変数として扱うカラム指定 categorical_features_indices = np.where(X.dtypes != np.float)[0]
訓練用と検証（Validation）用にデータを分割しておきます。
rom sklearn.model_selection import train_test_split X_train, X_validation, y_train, y_validation = train_test_split(X, y, test_size = 0.3)
モデリング
今回はクラス分類なので、catboost.CatBoostClassifierという関数を利用します。
import catboost mod = catboost.CatBoostClassifier(iterations=5000, calc_feature_importance = True, use_best_model=True, eval_metric = 'Accuracy' ) mod.fit(X_train, y_train, cat_features=categorical_features_indices, eval_set=(X_validation, y_validation), plot=True)
オプションについては、例えば以下のページにあります。
tech.yandex.com
今回は、CatBoostClassifierにはiterations：何回学習と検証を繰り返すか、calc_feature_importance：変数重要度を計算するか、use_best_model：iterationの中で一番検証用の指標が良いものを選ぶか、eval_metric：評価指標をオプションとして選んでいます。
モデル作成は、sklearnのようにfitメソッドで行います。cat_featuresでカテゴリカル変数として扱うカラムを指定し、eval_setで最初のvalidationに使用するデータを指定します。
実行後は、以下のような図がプロットされます。
これは、訓練・検証用データに対してAUCをプロットしたものです。点線が訓練用、普通の線が検証用なので、結構過学習しちゃってます（汗）。検証用データで一番良い場合はaccuracyが0.88程度のようです。
from sklearn.metrics import accuracy_score y_val_pred = mod.predict(X_validation) print('accuracy : %.2f' % accuracy_score(y_val_pred, y_validation)) # 出力 accuracy : 0.88
テストデータを用いて、予測してみます。これはsklearnの場合と同様に、predictメソッドで出力できます。ただし、floatとして出て来るので、出力後にintに変換します。こうしておかないと、Kaggleにsubmitしたときにscoreが0になります。
# テストデータはdf_raw_test df_raw_test.fillna(-999, inplace = True) y_test = mod.predict(df_raw_test) est = DataFrame() test['PassengerId'] = df_raw_test['PassengerId'] test['Survived'] = y_test # 予測値をintに変換 test['Survived'] = test['Survived'].astype('int') # 保存 test.to_csv('out/benchmark_catboost.csv', index = False)
この結果をKaggleにsubmitすると....
ほぼほぼ何も頑張ってないのにaccuracyが0.8を超えました。ただ、過学習しているモデルなので、より安定して精度を出したかったらパラメータ調整などをおこなって汎化性能を上げる必要があるかと思います。そういった努力は今後の課題にします。
変数重要度
Catboostはfeature_importances_で変数の重要度も出力することができます（ソートはされていない）。ざっくり可視化しておきます。
FI = DataFrame() FI = FI.append(mod.feature_importances_) FI.index = X.columns FI.columns = ['Feature_Importance'] # 可視化 import matplotlib.pyplot as plt %matplotlib inline FI.plot(kind='bar')
Pclassは乗客の社会的地位、Sexは性別、Ageは年代で、地位の高い女性や子供が優先的に救出されたという歴史的経緯が概ね反映されています。Ticletが効いているのはちょっと難しいですが。
Next Step
大雑把には使い方が分かったので、今後はGrid Searchなどを詰めていって、より使いこなせるようにしていこうと思います。
