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

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

Pythonでデータ分析:Prophetを使ってビットコインの予測(笑)をやってみる

導入

直近、これといって緊急の業務がなく、「自分の時間だ何勉強しようかなー」とPyStanとかをいじっていた矢先、「暇なら技術調査やってよ、Deep Learning的な何かとか」というお達しがきました。あいにく私は天邪鬼なので、2つ返事をして気になっていたけど触っていなかったProphetを調べることにしたのでした。

注:仕事はちゃんとしました(Seq2Seqの論文や書籍見て簡単な実装をしました)。

Prophet

Facebookが出した時系列予測のツールです。
facebook.github.io
すでに様々な方が紹介をしたり、Contributeしていたりするので、釈迦に説法感がありますが、このツールの良い点は、簡単に(分析の専門知識がなくても)ある程度それらしい予測値を出してくれるところです。ビジネス側でデータを活用したい場合や、分析者でもいったん簡単にデータから言えることを見てみる、といった場合に便利そうです。

というわけで、分析者としても抑えておきたいと前々から思っていたので、いい機会だと思って少し動かしてみました。一通りサンプルを見たあと、題材として何かあるかなと思って探していたところ、ビットコインのデータをPythonから取得できるライブラリを見つけたので、ビットコインを題材とすることにしました。

参考

  • Giihub(メソッドの引数とか見るのに使った)

github.com

  • ホクソエムさんの資料(最初に超見た)

d.hatena.ne.jp

darden.hatenablog.com

データ取得

以下のコードを使います。
github.com
インストール方法はサイトに書かれています。pipが使えるなら

$ pip install https://github.com/s4w3d0ff/python-poloniex/archive/v0.4.6.zip

とすれば使用できる状態になるはずです。

データですが、以下のコードを実行すると、ビットコイン(US)のデータ(日次)を過去500日分取得してくれます。ざくっと紹介すると、polo.returnChartDataでデータを辞書型で取得しています。periodでどれくらいの間隔のデータかを指定(5分足なども指定できます)しています。ただ、時間がUNIX timeからの経過時間(秒)になっているので、後でよしなに変換しています。

# numpyやpandas
import numpy as np
import pandas as pd
from pandas import DataFrame, Series

# APIに関係するライブラリ
import poloniex
import time
import datetime

def getDataPoloniex():
    polo = poloniex.Poloniex()
    polo.timeout = 2
    chartUSDT_BTC = polo.returnChartData('USDT_BTC', period=polo.DAY, start=time.time() - polo.DAY * 500, end=time.time())
    tmpDate = [chartUSDT_BTC[i]['date'] for i in range(len(chartUSDT_BTC))]
    date = [datetime.datetime.fromtimestamp(tmpDate[i]).date() for i in range(len(tmpDate))]
    data = [float(chartUSDT_BTC[i]['open']) for i in range(len(chartUSDT_BTC))]
    return [date, data]

# データ取得
dat = getDataPoloniex()
# DataFrameへ格納
df = DataFrame(dat).T
df.columns = ['ds', 'y']

ちなみに、プロットすると以下のようになっています。

import matplotlib.pyplot as plt
%matplotlib inline
df.plot(x = df.columns[0], figsize = (15, 8))

f:id:tekenuko:20171018000038p:plain
途中で下落したりすれど、急成長しているような系列になっています。

Prophetによるモデル作成

ProphetはPythonの場合はpipで導入することができます。

pip install fbprophet

PyStanとか色々依存ライブラリがあるようなんですが、そういったものたちは先に入れとくといいのかもしれません。自分はすべて先に入っていたので比較ができませんが。

モデル構築はsklearnライクにできます。

mod = Prophet()
mod.fit(df)

非線形トレンドを入れたり、変化点を考慮したりといろいろオプションはあるのですが、ここではよしなにできる感を演出するのに何も指定しないでおきます。線形トレンド以上の成長をしていそうな系列なので、ホントはちゃんと処理したいですけどね。

予測に関しては、make_future_dataframeとpredictで新たに生成した系列(100日分)に予測値を格納します。

future_df = mod.make_future_dataframe(periods = 100)
forecast_df = mod.predict(future_df)

簡単なプロットはProphetじたいにplotというメソッドがあり、それでよしなな図が見れます。

mod.plot(forecast_df)
plt.show()

f:id:tekenuko:20171018000907p:plain

まあ、線形トレンド的な予測がされています。もっと長期まで予測しようと思ったら大まかな傾向(指数的)とは異なってくるとは思いますが、おおかなに上昇していきそうというのが捉えられている、としましょう。

もう少し細かい成分でプロットする機能もあります。

mod.plot_components(forecast_df)

f:id:tekenuko:20171018001512p:plain

週の上下している感を見ると

  • 火曜日は上がり調子
  • 水曜、木曜にかけて下がる

というホントかな、、という傾向が見られます。ちょっとこれは期間を変えて検証を繰り返す、とかしないと正確なことは何も言えません。ちなみに、他の記事では学習期間を変えると週の振る舞いは変わると言っているものもありますので、やはりこれから何かを主張するのは難しいのかもしれません。

qiita.com

感想

データがあるときに、チャチャッと分析をするのに便利です。これからは割りと使用するかもしれない、そう思われるツールでした。ただ、やはり適用限界はありそうなので、100%頼るというよりはできそうなところまで使ってみる、という付き合い方になるかなと思います。

Pythonでデータ分析:主成分分析(PCA)による異常検知

導入

データ分析の種類の一つとして、教師なし学習による異常検知というものがあります。ほとんどが正常なデータでまれに異常なデータが混じっている、その異常発生のパターンや異常と他の要因との紐付きがいまいちつかみきれていないというような場合、教師あり学習による2値分類がうまくワークしない、といった状況がありえます。そういった場合には、正常パターンを教師なし学習で学び、その正常パターンから外れているものを異常とする、という方法が有効です。

異常検知の方法ですが、正常パターンを正規分布でフィットするようなものから状態空間モデルに即したものまで多くのバリエーションがあります。その中の一つに、主成分分析(PCA)をベースにした異常検知の方法があります。正常データがある程度相関を持っているような多変数正規分布に従っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています(正規分布から外れている場合は、Box.Cox変換などで無理やり正規分布に近づける、他の手法を使うなどを検討します)。この方法ですが、Rでは実行例は井出さんの異常検知の本
www.coronasha.co.jp
にありますが、Pythonでの実行例というのが書籍やネットを探してもいまいち見つかりません。そのため、メモがてらこちらに簡単に残しておこうと思います。

多変数正規分布による人工データ

異常検知のデモ用に、人工データを用意します。PCAによる異常検知の場合、データは正規分布に従い、かついくつかの変数で相関が高いような場合にうまくはまります。というわけで、多変数正規分布で変数間にいくらか相関のあるようなデータを生成します。

# 必要なライブラリ
import numpy as np
import pandas as pd
from pandas import DataFrame, Series

# 多変数正規分布に従う乱数を生成
mean = (4, 2, 3)
sigma = [[1, 0.7, 0.75], [0.7, 1, 0.3], [0.75, 0.3, 1]]
df = np.random.multivariate_normal(mean, sigma, 500)

主成分分析(PCA)

※以下のコードは、無駄にDataFrameが使われていますが、もっとすっきり書ける可能性があります。

主成分分析は、次元圧縮の方法の一つです。データセットに含まれる変数が非常に多いがいくらか相関が強いものがあるとき、適切に基底の回転・取り替えを行うことによってより少ない軸でデータを近似することが期待できます。主成分分析(PCA)は今のデータセットから新しい基底(主成分)を求め、その基底のうちデータの情報を失わない(軸に射影しても分散が小さくならない)ものでいくらか代表して近似してしまおう、といった手法になっています。

PCAは、sklearnのメソッドで簡単に行うことができます。

# 異常度算出の都合上、正規化を行う
std = StandardScaler()
X = DataFrame(std.fit_transform(df))

from sklearn.decomposition import PCA
# インスタンス化
pca = PCA()
# 最初の1/5を学習に用いる
pca.fit(X[0:round(X.shape[0] / 5)])

pca.fitでレコード数×特徴量数の行列に関して特異値分解(もしくはこの行列と転置をかけた共分散行列というものの固有値分解)を行い、新しい基底を求めています。今回、変換前の軸は3つあったのですが、変換後は2つの軸でデータを近似することにします(こちらに関しては、本来は軸の寄与率などでどれくらい軸を採用するかを決めるなどします)。新しい基底はcomponents_というメソッドで引き出すことができ、その行方向を制限します。

Um = DataFrame(pca.components_[0:2, :]).T

ここで、後の都合上転置をしています。また、共分散行列というもの( S^2 = {I}_2 - U_m U_m^T)を以下で定義しておきます。

S2 = np.identity(X.shape[1]) - Um.dot(Um.T)

異常度算出

PCAで次元圧縮した際に用いたデータは、正常なデータのみ、もしくはほとんど異常なデータは存在しないとします。すると、正常データは近似した軸で貼られる部分空間上でうまく近似でき、逆に異常なデータはその部分空間から大きく離れるのではないかと期待されます。この部分空間は正常部分空間と呼ばれます。異常度を評価する一つの方法が、データと正常部分空間からの距離を用いる方法です。今、正規化(中心化でもよい)されたデータを\vec{x}(次元は特徴量の数)とすると、異常度a(\vec{x})は次のように書かれます。
$$a(\vec{x}) = \vec{x}^T S^2 \vec{x}^T$$
導出は井出さんの異常検知の本などを参照してください。


Pythonでは、おそらくライブラリで異常度を出力してくれるメソッドはないので、手で書く必要があります。ここでは、1レコードに対して異常度を計算する関数を作っておき、それを検証用データに対してループさせて使う、としています。

1レコードに対して異常度算出
def CalcAnomaly(x, S2):
    return x.T.dot(S2).dot(x)

# 検証データ(全データ)に対して異常度算出
a = [CalcAnomaly(X.iloc[i, :], S2) for i in range(X.shape[0])]   

検証・可視化

各レコードに対し、異常度を可視化してみます。

import matplotlib.pyplot as plt
%matplotlib inline

plt.figure(figsize = (15, 8))
plt.plot(a)
plt.xlabel('n_step')
plt.ylabel('anomaly score')

f:id:tekenuko:20171016205305p:plain

概ね1より小さい値になっています。まれに1程度の相対的に異常度が高い値が出力されていますが、こういった値は乱数の関係上発生しうるものだと思われます。本当に異常なデータがある場合は、今のスケール感でいうと数十といった値が出てくることもあるので、今回はほぼ正常データしかないだろうと思われるデータになっています。

また、もとの空間での散布図と異常度の関係も可視化しておきましょう。3軸なので、散布図も3次元でプロットしています。

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize = (12, 10))
ax = Axes3D(fig)
p = ax.scatter(df[:, 0], df[:, 1], df[:, 2], marker = 'o', c = a, cmap = 'jet', s = 50)
fig.colorbar(p)

f:id:tekenuko:20171016205724p:plain

新しい軸はデータの相関のある方向(図で言うと斜め方向)などになっているのですが、その軸たちから離れているデータほど異常度が大きくなっている、といった様子が図から読み取れます。このように、PCAによる異常検知は正常データがある程度相関を持っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています。

Next Step

今回は、主成分分析(PCA)をベースにした異常検知の方法があります。PCAによる異常検知は正常データがある程度相関を持っており、異常データはその相関関係が崩れたようなものになっている、といった場合に有効な方法になっています。異常検知の方法は他にも色々あるので、機会があればまた別の方法も紹介したいと思います。

Memo:Gluonの解説やコード紹介(海外)

導入

2017年10月12日(現地時間)に、MicrosoftAWSがGluonというDeep Learningのライブラリを公開しました。

www.itmedia.co.jp

日本語の解説記事があまり見当たらなかったので、簡単なところは自分で試してみるなどし、いくつか記事にもしました。色々調べていたところ、Documentや国際会議などではすでにGluonじたいやコードの紹介があることがわかってきました。こういった情報を整理しておくのは後々の自分にとって有用だろうと思ったので、メモがてら残しておきます。

2017年の国際会議でのTalk

CVPR2017

mli.github.io

ICML2017, KDD2017

github.com


GithubのLecture Note

ICML2017, KDD2017での発表者のGithubにMXNetも含めたLecture Noteがあります。形式はJupyter Notebookです。MXNetのメソッドの使い方から、深層強化学習まで、広い範囲を扱っています。これを最初に見つけてたらどんなに楽だったか。
github.com

また、上記のLectureに対応するWeb Pageもあります。

所感

日本では、ついこの間登場した感が出てますが、海外の情報や国際会議のTalkではPreliminalなものも含め前々から紹介されていたりする、そんな当たり前のことを再確認しました。日本人だと日本語の記事を情報源にしたくなりますが、真面目に情報収集したいならソース元を探しにいかないと。論文のサーベイと同じですね。

とはいえGluonは海外も含めてもまだまだ情報が少ない状況なので、情報収集は定期的に続けていき、それを自分が紹介することで日本での紹介や実例がより活発になっていけばよいなと思います。

GluonでDeep Learning:CNNを組んでみる

導入

前回、MicrosoftAWSが公開したライブラリである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日(現地時間)に、MicrosoftAWSが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

インストール

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でロジスティック回帰なんかが具体的にやろうと思っていることです。

Pythonでデータ分析:Catboost

導入

2017年7月に、ロシアのGoogleと言われている(らしい)Yandex社から、Catboostと呼ばれるGradient Boostingの機械学習ライブラリが公開されています。

catboost.yandex

ここ何ヶ月か調整さんになっていて分析から遠ざかりがちになりやすくなっていたのですが、手を動かしたい欲求が我慢できなくなってきたので、いい機会だと思って触ってみることにしました。

インストール

anacondaなどでPythonを導入していることを前提にします。その場合、おそらくpipで一瞬でインストールできます。

$ pip install catboost

使用データセット: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に使用するデータを指定します。

実行後は、以下のような図がプロットされます。
f:id:tekenuko:20171013225210p:plain
これは、訓練・検証用データに対して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すると....

f:id:tekenuko:20171013230117p:plain

ほぼほぼ何も頑張ってないのに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')

f:id:tekenuko:20171013230525p:plain

Pclassは乗客の社会的地位、Sexは性別、Ageは年代で、地位の高い女性や子供が優先的に救出されたという歴史的経緯が概ね反映されています。Ticletが効いているのはちょっと難しいですが。

Next Step

大雑把には使い方が分かったので、今後はGrid Searchなどを詰めていって、より使いこなせるようにしていこうと思います。