プロクラシスト

今日の寄り道 明日の近道

【Day-20】PyTorchを超絶使いやすくするsklearnラッパー『skorch』で快適ディープラーニング

f:id:imslotter:20171220124248p:plain

データ分析ガチ勉強アドベントカレンダー 20日目。

Skorchとは

f:id:imslotter:20171220123959p:plain
  • PyTorchのsklearnラッパー
  • sklearnのインターフェースでPyTorchを使える

自分が感じたメリット

  • sklearnが行うところとPyTorchが作るところがはっきりしていて、コードがすっきりまとまる
  • sklearnの関数(fit, predictなど)がそのまま使える
  • numpy配列をそのまま学習器に投げられるので、イメージしやすい
  • 前処理、GridSearch, 評価もsklearnのものを使えて、便利
  • インストール

    git clone https://github.com/dnouri/skorch.git
    cd skorch
    # create and activate a virtual environment
    pip install -r requirements.txt
    # install pytorch version for your system (see below)
    python setup.py install
    

    ※私の環境(ubuntu14.04, python3.6.1)では、sklearn0.18.1では実行時エラーが出たので、sklearn0.19.1にアップグレードしました。

    pip install -U sklearn

    使い方

    データ読み込みはsklearn

    データの作り方などはday-7の記事を参照

    コードをみる
    %matplotlib inline
    import matplotlib.pyplot as plt
    
    import torch
    from torch import nn
    import torch.nn.functional as F
    torch.manual_seed(0);
    
    import numpy as np
    from sklearn.datasets import make_classification
    import sklearn
    X, y = make_classification(1000, 20, n_informative=10, random_state=0)
    X = X.astype(np.float32)
    X.shape, y.shape,y.mean()
    

    学習ネットワークの構築はPyTorch

    day19の記事のように、ネットワークの設計をdefine by runでしていきましょう。

    コードをみる
    class ClassifierModule(nn.Module):
        def __init__(
                self,
                num_units=10,
                nonlin=F.relu,
                dropout=0.5,
        ):
            super(ClassifierModule, self).__init__()
            self.num_units = num_units
            self.nonlin = nonlin
            self.dropout = dropout
    
            self.dense0 = nn.Linear(20, num_units)
            self.nonlin = nonlin
            self.dropout = nn.Dropout(dropout)
            self.dense1 = nn.Linear(num_units, 10)
            self.output = nn.Linear(10, 2)
    
        def forward(self, X, **kwargs):
            X = self.nonlin(self.dense0(X))
            X = self.dropout(X)
            X = F.relu(self.dense1(X))
            X = F.softmax(self.output(X), dim=-1)
            return X
    

    skorchでwrap

    skorch.net 用意されているのは2種類

    関数 用途
    NeuralNetClassifier 分類器をsklearn風に
    NeuralNetRegressor 回帰をsklearn風に

    初期化の際に、学習の仕方を決める。パラメータは、PyTorchの関数を使用できる。

    • criterion : 損失関数の設定
    • optimizer : 最適化関数の設定
    • lr : 学習率の決定

    その他色々とパラメータはあるので、skorch.netを参照。

    PyTorchのパラメータに関しては『PyTorch入門』 使い方&Tensorflow, Keras等との違いとは?を参照

    • .fitで、自動的にtorch.tensorに変換される
    • validationまでやってくれて嬉しい
    コードをみる
    from skorch.net import NeuralNetClassifier
    net = NeuralNetClassifier(
        module=ClassifierModule,
        max_epochs=20,
        lr=0.1,
        # use_cuda=True,  # uncomment this to train with CUDA
    )
    net.fit(X,y)
    y_pred = net.predict(X[:5])
    y_proba = net.predict_proba(X[:5])
    for pred, proba in zip(y_pred,y_proba):
        print("score {}: class {}".format(proba, pred))
    

    sklearnとのその他連携

    pipeline

    • スケーリングなどsklearnの処理をデータの流れに組み込める
    • 前処理の話は day-8
    コードをみる
    from sklearn.pipeline import Pipeline
    from sklearn.preprocessing import StandardScaler
    
    # データの流れるPipelineを設計
    pipe = Pipeline([
        ("scale",StandardScaler()),
        ("neuralnet", net)
    ])
    print(pipe.named_steps)
    
    pipe.fit(X,y)
    
    y_pred = pipe.predict(X[:5])
    y_proba = pipe.predict_proba(X[:5])
    for pred, proba in zip(y_pred,y_proba):
        print("score {}: class {}".format(proba, pred))
    
    • NeuralNetClassifier/NeuralNetRegressorで設定できるパラメータを調べることができる(GridSearchCV or RandomSearchCV)
    • optimizerなどもtorch.optimのリストを使って網羅的に調べられる
    • クラス内変数も module__hoge(hogeはメンバ変数)と入力することで置き換えられる
    コードをみる
    from sklearn.model_selection import GridSearchCV
    params = {
        "lr":[i*0.01 for i in range(1,5)],
        "optimizer":[torch.optim.Adam, torch.optim.Adagrad, torch.optim.SGD],
        "module__num_units":[10,20],
    }
    gs = GridSearchCV(net, params, refit=False, cv=3, scoring="accuracy")
    gs.fit(X,y)
    import pandas as pd
    df = pd.DataFrame(gs.cv_results_)
    df_scored = df.sort_values(by=["rank_test_score"])[["params","mean_test_score","std_test_score","mean_fit_time"]]
    df_scored
    

    MNIST

    MNISTを実装してみる。sklearnで書けるところはsklearnで書いて、PyTorchで書くべきところはそちらで書く。使い分けがはっきりしていてかなりいい感じ。

    # データの読み込み(sklearn)
    from skorch import NeuralNetClassifier
    from torch import nn
    from sklearn.datasets import fetch_mldata
    from sklearn.model_selection import train_test_split
    from sklearn.metrics import classification_report
    mnist = fetch_mldata('MNIST original')
    X = mnist.data.astype('float32')
    y = mnist.target.astype('int64')
    X /= 255
    XCnn = X.reshape(-1, 1, 28, 28)
    XCnn_train, XCnn_test, y_train, y_test = train_test_split(XCnn, y, test_size=0.25, random_state=42)
    # Networkの設計(PyTorch)
    class Net(nn.Module):
        def __init__(self):
            super(Net, self).__init__()
            self.conv1 = nn.Conv2d(1, 32, kernel_size=3)
            self.conv2 = nn.Conv2d(32, 64, kernel_size=3)
            self.conv2_drop = nn.Dropout2d()
            self.fc1 = nn.Linear(1600, 128) # 1600 = number channels * width * height
            self.fc2 = nn.Linear(128, 10)
    
        def forward(self, x):
            x = F.relu(F.max_pool2d(self.conv1(x), 2))
            x = F.relu(F.max_pool2d(self.conv2_drop(self.conv2(x)), 2))
            x = x.view(-1, x.size(1) * x.size(2) * x.size(3)) # flatten over channel, height and width = 1600
            x = F.relu(self.fc1(x))
            x = F.dropout(x, training=self.training)
            x = self.fc2(x)
            x = F.softmax(x, dim=-1)
            return x
    # ラッパーを使う(skorch)
    net = NeuralNetClassifier(
        Net,
        max_epochs=10,
        lr=1,
        optimizer=torch.optim.Adadelta,
        # use_cuda=True,  # uncomment this to train with CUDA
    )
    # training
    net.fit(XCnn_train, y_train)
    
    # test
    y_pred = net.predict(XCnn_test)
    print(classification_report(y_test, y_pred))
    

    結果

    学習の様子

    結果の表示

                 precision    recall  f1-score   support
    
              0       1.00      0.99      0.99      1677
              1       0.99      0.99      0.99      1935
              2       0.99      0.99      0.99      1767
              3       0.99      0.99      0.99      1766
              4       0.99      0.99      0.99      1691
              5       0.99      0.99      0.99      1653
              6       0.99      0.99      0.99      1754
              7       0.99      0.99      0.99      1846
              8       0.98      0.99      0.99      1702
              9       0.99      0.98      0.99      1709
    
    avg / total       0.99      0.99      0.99     17500
    
    

    まとめ

    sklearnとPyTorchをつなげられる便利なライブラリを紹介しました。 PyTorchの持つ柔軟性と、sklearnの利用しやすさが相まって最高ですね。

    かなり使いやすいので、sklearnに慣れている人は、一回使ってみることをおすすめします!

    明日はどうするかな、ベイズ的な変化点検知をする予定だが、果たして、、、ではでは!

    【Day-19】『PyTorch入門』 使い方&Tensorflow, Keras等との違いとは?

    f:id:imslotter:20171219000627p:plain

    データ分析ガチ勉強アドベントカレンダー 19日目。

    2日間、Kerasに触れてみましたが、最近はPyTorchがディープラーニング系ライブラリでは良いという話も聞きます。

    とりあえずTutorialを触りながら使ってみて、自分が疑問に思ったことをまとめていくスタイルにします。

    また、同じく有名ライブラリであるKerasやTensorFlowについての比較もしたいと思っています(Define and RunかDefine by Runか)

    PyTorchとは

    f:id:imslotter:20171219043523p:plain

    ざっくりまとめると

    • 深層学習用ライブラリ
    • chainerからフォークされて作られている。思想はChainerを引き継いでいる。
    • GPUでの計算も楽に行える
    • Define by Run
    • コミュニティが活発で、後発にもかかわらず急速に発達

    特にDefine by RunはTensorflowやKerasと大きく異なる特徴です(後述)

    PyTorch入門

    変数の扱い方

    PyTorchでは、変数や関数をPyTorchオリジナルの型として扱う。大体はnumpyライクに扱えるので、ストレスは少ない

    x = torch.Tensor(5,3) #ゼロ行列
    y = torch.rand(5,3) # ランダム行列
    # 演算
    print(x.add_(y)) #アンダースコアをつけると破壊的変化
    print(torch.add(x,y)) #torchの関数から
    print(x+y) #通常の演算もいける
    # numpy likeに扱える
    print(y[:,1],y.mean(), y.std())  #ある程度の関数はそろっている
    print(torch.ones(10))
    print(torch.ones(10)+1)
    # numpyからの変換も可能
    import numpy as np
    a = np.arange(4).reshape(2,2)
    x = torch.from_numpy(a)
    

    Autograd

    • torch.autograd FunctionVariablesで定義すると、演算結果が保存されるので、 そのままさかのぼった微分が可能。

    下記の式は、

    • {x=} [[1,2],[3,4]]
    • {y = x^{2} +2}
    • {z = y^{2} = (x^{2}+2)^{2}}

    {out = \frac{1}{4} \sum_{i} z_{i} z_{i}  (x^{2}+2)^{2}} より、 {\frac{\partial out}{\partial x_{i}}=x_{i}(x_{i}^{2}+2)}

    その結果を返す より詳しい内容は公式

    from torch.autograd import Variable
    import numpy as np
    x = Variable(torch.FloatTensor([1,2,3,4]), requires_grad=True)
    y = x.float()**2 + 2
    z = y**2
    out = z.mean()
    out.backward() # 
    print(x.grad)
    

    出力

    <AddBackward0 object at 0x7f88d576ea90>
    Variable containing:
      3  12
     33  72
    [torch.FloatTensor of size 2x2]
    

    チュートリアル:NeuralNetの構築

    学習の手順

    • 学習ができるニューラルネットを定義する。
    • データセットを繰り返し投げられるようにする。
    • インプットに応じて、処理(計算)が走るようにする
    • 誤差(損失)を計算する
    • back propagationを実施する
    • optimizerによってパラメータを更新する

    この中でデータセットを繰り返し投げる以外の要素を説明していく。

    ライブラリ構成

    PyTorchのよく使うパッケージ

    パッケージ 説明
    torch NumPyのような配列(テンソル)ライブラリ。GPUも簡単に使える
    torch.autograd 自動微分ライブラリ。Torchで定義されたすべての関数に対して微分を可能にしている
    torch.nn ネットワークを構成する要素
    torch.optim パラメータ更新時に使う最適化パッケージ
    torch.utils DataLoaderなど、その他のユーティリティ関数

    torch.nn

    ネットワークを構成する要素関連はここに入っている

    • Conv1d,Conv2d,MaxPool1d, AvePool1d, AdaptiveAvgPool1d... : CNN向け
    • RNN, LSTM, GRU.... : RNN向け
    • BatchNorm1d,BachNorm2d, InstanceNorm1d,... : バッチ正規化
    • functional.hoge : 活性化関数(SoftMax, ReLUなど)
    • その他、dropout, sparse_embeddings,...

    詳しくはpytorch:nn

    すべてPyTorch内の関数で渡してやることができれば、autogradにより、back propagation時の微分演算などがスムーズ

    import torch
    from torch.autograd import Variable
    import torch.nn as nn
    import torch.nn.functional as F
    
    class Net(nn.Module):
    
        def __init__(self):
            super(Net, self).__init__()
            # 1 input image channel, 6 output channels, 5x5 square convolution
            # kernel
            self.conv1 = nn.Conv2d(1, 6, 5)
            self.conv2 = nn.Conv2d(6, 16, 5)
            # an affine operation: y = Wx + b
            self.fc1 = nn.Linear(16 * 5 * 5, 120)
            self.fc2 = nn.Linear(120, 84)
            self.fc3 = nn.Linear(84, 10)
    
        def forward(self, x):
            # Max pooling over a (2, 2) window
            x = F.max_pool2d(F.relu(self.conv1(x)), (2, 2))
            # If the size is a square you can only specify a single number
            x = F.max_pool2d(F.relu(self.conv2(x)), 2)
            x = x.view(-1, self.num_flat_features(x))
            x = F.relu(self.fc1(x))
            x = F.relu(self.fc2(x))
            x = self.fc3(x)
            return x
    
        def num_flat_features(self, x):
            size = x.size()[1:]  # all dimensions except the batch dimension
            num_features = 1
            for s in size:
                num_features *= s
            return num_features
    
    net = Net()
    print(net)
    

    出力 : ネットワーク構成が分かる

    Net(
      (conv1): Conv2d (1, 6, kernel_size=(5, 5), stride=(1, 1))
      (conv2): Conv2d (6, 16, kernel_size=(5, 5), stride=(1, 1))
      (fc1): Linear(in_features=400, out_features=120)
      (fc2): Linear(in_features=120, out_features=84)
      (fc3): Linear(in_features=84, out_features=10)
    )

    値を投げる

    torchの型にしたがってデータを構成し、投げる。

    params = list(net.parameters())
    input = Variable(torch.randn(1, 1, 32, 32))
    out = net(input)
    

    損失関数

    誤差の算出のために、いろいろな関数が用意されている。シンプルなのはMean Squared Error。いろんな誤差を知りたい人は、機械学習で使う指標総まとめ(教師あり学習編)へ。

    • MSELoss
    • CrossEntropyLoss
    • NLLLoss(negative log likelihood loss) (1d, 2d)
    • Poisson NLL Loss
    • KLDivLoss (Kullback-Leibler divergence)
    • BCELoss (Binary Cross Entropy)
    • ...

    より詳しくは公式

    output = net(input)
    target = Variable(torch.arange(1, 11))  # a dummy target, for example
    criterion = nn.MSELoss()
    loss = criterion(output, target)
    

    back propagation

    損失が算出できたので。この損失からさかのぼっていく。PyTorchだと、backwardとだけ実行すれば、簡単にパラメータの微分値を算出できるようになっている。

    net.zero_grad() #とりあえず0で初期化
    loss.backward()
    

    重みのアップデート

    学習には、torch.optimを使う。下記のようなアルゴリズムがある。

    • SGD, ASGD, Nesterov-SGD
    • Adadelta, Adagrad, Adam, SparseAdam, Adamax
    • LBFGS
    • RMSProp, Rprop

    また、torch.optim.lr_scheduler学習率の調整もできる。

    import torch.optim as optim
    
    # create your optimizer
    optimizer = optim.SGD(net.parameters(), lr=0.01)
    
    # in your training loop:
    optimizer.zero_grad()   # zero the gradient buffers
    output = net(input)
    loss = criterion(output, target)
    loss.backward()
    optimizer.step()    # Does the update
    

    これで一通り終了。最終的にはデータを繰り返し投げて、学習を進めていく。 Githubのday19では、最後にMNISTの分類をPyTorchでやっている。そちらも是非どうぞ。

    主要フレームワークとの比較とDefine by Run

    Define and RunとDefine by Run

    ディープラーニングライブラリは上記の二種類の設計思想に大きく分かれる。 Chainerの血を受け継ぐPyTorchは、Define by Run。TensorflowやKerasとの最も大きな違い。

    どのように違うのか。下図は、PyTorchとKerasで定義した、Mnistに投げるCNN。 PyTorchがデータを投げて実行しながらネットワークを定義するのに対して、Kerasはネットワーク構成をかっちりと決めてからデータを投げる。定義の時点でデータは考えない。

    Define by Runのメリットとしては、

    • ネットワークを動的に変更できる。

    ということで、研究者界隈では、Define by Runが主流。

    使いやすさのKeras, 柔軟性のPyTorch

    とはいえ、Kerasはやはり圧倒的に使いやすい。 ネットワークを積み上げるだけで、簡単にディープラーニングができてしまう。

    その一方で、ネットワークの詳細設計はしにくい。

    通常のデータ分析業務ならば、Kerasの簡潔さが輝くだろうし、研究や難しいタスクならばPyTorchが優位なのかなという印象。

    コミュニティの強いPyTorch

    PyTorchの優位性は、コミュニティの強さにもある。 同じくdefine by runのChainerは、やはり日本が主流であって、海外での利用は少ない。

    一方でPyTorchは、研究者が続々と採用していて、主要な論文が出るとすごいスピードで実装されるので、最新論文のアルゴリズムをすぐに試すことができる。the incredible pytorchなどが好例。

    github.com

    まとめ

    PyTorchについて勉強してみた。 設計思想の違いなんかにも触れられたと思う。一長一短なので、両方とも使えるようになっていたら便利そうだなと思った。

    明日もPyTorch系の話題。どうやらsklearnラッパーがあるようなので、試してみたいと思う。ではでは!

    詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

    詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

    Natural Language Processing With Pytorch: Build Intelligent Language Applications Using Deep Learning

    Natural Language Processing With Pytorch: Build Intelligent Language Applications Using Deep Learning

    【Day-18】時系列のディープラーニング、RNNのまとめとKeras実装

    f:id:imslotter:20171218194346p:plain

    データ分析ガチ勉強アドベントカレンダー 18日目。

    Kerasの使い方を復習したところで、今回は時系列データを取り扱ってみようと思います。 時系列を取り扱うのにもディープラーニングは用いられていて、RNN(Recurrent Neural Net)が主流。

    今回は、RNNについて書いた後、Kerasで実際にRNNを実装してみます。

    今日の参考書はこの本!

    詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

    詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

    • 作者:巣籠 悠輔
    • 発売日: 2017/05/30
    • メディア: 単行本(ソフトカバー)

    RNNとは

    通常のNeural Networkとの違い

    Recurrent Neural Netの略。再帰ニューラルネット。何が再帰かは下図

    入力層、隠れ層、出力層の構造はDay-15で説明したNeural Netとほとんど同じだが、最も大きな違いは隠れ層同士がつながっているところ。これによって時系列に対応できるようにしている。順伝播の計算も、再帰性の要素が入ってくる。(添え字等はDay15参照)

    誤差逆伝播法のアルゴリズム

    基本的にニューラルネットと同じだが、隠れ層の時間方向に微分がつながっていくので、時間を遡りながら誤差の影響を伝えていく、しかし、本来は一番初めの時刻まで遡るものだが、遡りすぎると時間方向の微分がかけ粟さて勾配消失/爆発を引き起こし、うまく学習されないことも多くある。シンプルなRNNではBPTTを途中で打ち切る場合も多い。

    下記にアルゴリズムを記す

    1. 順伝播計算により{z}の更新式、出力{y}を求める

    2 ラベル{d_{nk}^{t}}との誤差を計算

    3. 重みを最小化する


    4. デルタはBPTTを行うことで再帰的に求められる

    5. 何らかのオプティマイザにより重みを更新

    勾配消失の工夫 : LSTMやGRU

    時間を遡ってBackpropagationを行うため、勾配消失が起こる。それを解消する工夫として、LSTM(Long-Short Term Memory)とGRU(Gated Recurrent Unit)がある。上の図の隠れ層に置ける、○の中身を工夫しながら、何を忘れて何を覚えておくか**まで賢く学習する。

    LSTM

    正直、本の図はちょっと分かりにくい、ウェブをあさってて一番厳密かつ分かりやすかったのは A Beginner’s Guide to Recurrent Networks and LSTMsというサイトの図

    複雑...一番のポイントは

  • 中心のcellと書かれている箇所(CECという)+forget gateで、過去情報をどのくらい記憶するかを調整する
  • これにより、BTPP中の勾配消失を数式的にも抑えることができた*1

    GRU

    LSTMの難点は複雑でパラメータが多いところ。上図で言えばinput/forget/blockと、それぞれのGateに入れるときの重みを最適化しないといけない。一方でGRUはLSTMベースだが、入力ゲートがLSTMよりも断然少ない。

    ゆえに計算量も少ない。さらに、タスクによってはLSTMより良い性能を発揮するということで、注目されている

    見た目の違いが論文に載っている。(Empirical Evaluation of Gated Recurrent Neural Networks on Sequence Modeling)

    だいぶシンプル。計算も軽いわけだから、とりあえずGRUから試してみるのがいいのかもしれない。

    自然言語をにぎわすAttention Model

    NLP(Natural Language Process)のRNN界隈では、Attention Modelというのがにぎわっているように見受けられる。

    勉強不足でまだまだ全然分かっていないので、勉強用のソースだけ張っておく。勉強。

    Keras実装

    練習がてら、KerasでRNNを実装してみる。 githubにもあげているのでそちらも参考に。

    github.com

    データの入力の仕方がちょっと難しかったので、そこだけ重点的に。

    データの変形、入力

    BPTTの都合もあり、投げるデータの配列をすこしいじってやらないといけない。いじり方について図でまとめてみた。このような3次のテンソル的な配列になる。

    コードをみる
    %matplotlib inline
    import matplotlib.pyplot as plt
    import pandas as pd
    import numpy as np
    import keras
    df = pd.read_csv("AirportPassengers.csv",delimiter=";").dropna()
    data = []
    target = []
    max_len = 24
    dim = 1
    # 正規化
    maximum = df.Passengers.max()
    minimum = df.Passengers.min()
    df["Passengers"] = (df.Passengers-minimum)/(maximum-minimum)
    # データを箱に入れる
    for i in range(len(df)-max_len-1):
        data.append(df.Passengers.values[i:i+max_len])
        target.append(df.Passengers.values[i+max_len+1])
    # データの整形
    data = np.array(data).reshape(len(data),max_len,dim)
    target = np.array(target).reshape(-1,1)
    # データの分割
    from sklearn.model_selection import train_test_split
    N_train = int(len(data)*0.7)
    N_test = len(data) - N_train
    X_train, X_validation, Y_train, Y_validation = train_test_split(data, target, test_size=N_test)
    

    学習モデル構築

    keras.layers.recurrentの中に、学習モデルが入っている。 LSTM, GRUも実装されていて、model.addするだけで使えるようになる。Optimizerの設定は普通のニューラルネットと同じ。Day-17にOptimizerについて書かれている

    コードをみる
    from keras.models import Sequential
    from keras.layers import  Dense
    #一番シンプルなやつ
    from keras.layers.recurrent import SimpleRNN
    from keras.layers.recurrent import LSTM, GRU #改良版
    model = Sequential()
    # ネットワーク構成
    input_shape=(max_len, dim)
    # model.add(SimpleRNN(units=64, kernel_initializer="random_uniform",input_shape=input_shape))
    #model.add(LSTM(units=64, input_shape=input_shape))
    model.add(GRU(units=64, input_shape=input_shape))
    
    model.add(Dense(input_shape[1],activation="linear"))
    # optimizerの設定
    from keras.optimizers import Adam
    model.compile(loss="mse", optimizer=Adam())
    

    学習と予測

    7割を学習用データにして、残りの3割を予測することにした。

    コードをみる
    from keras.callbacks import EarlyStopping
    early_stopping = EarlyStopping(monitor='val_loss', patience=10, verbose=1)
    epochs = 500
    batch_size = 10
    model.fit(X_train, Y_train,
              batch_size=batch_size, epochs=epochs,
              validation_data=(X_validation, Y_validation))
    prediction = []
    data_in = data.reshape(data.shape[0],data.shape[1])[N_train]
    
    iteration = len(data) - N_train
    for _ in range(iteration):
    #     print(data_in)
        pred = model.predict(data_in.reshape(1,-1,1))
        data_in = np.delete(data_in, 0)
        data_in = np.hstack((data_in, pred[0]))
        prediction.append(pred[0,0])
    passenger = list(df.Passengers)
    data_num = len(passenger)
    plt.plot(passenger)
    plt.plot(range(N_train+max_len, N_train+max_len+iteration), prediction)
    

    結果

    上記の学習を3回行ってみて、LSTMと普通のRNNとGRUを比べた。

    今回の結果は、LSTMが良かった。

    きちんとトレンド/周期共に捕らえられているようだった。ただ、時間を追うごとに値がずれていっている。もう少し精度を上げたいなら、前処理をしっかりしたほうがいいのかもしれない。前処理については以前時系列まとめで書いた(【Day-12】時系列分析の良リソースまとめ&基礎チュートリアル)

    まとめ

    せっかくKerasを学んだので、時系列もディープでしておこうと思い、まとめた。 ちょっと苦労したけど、これで学習器に入れられるようになった。予測がちゃんとできれば楽しい。

    明日はPyTorchやりたいなぁ。ではでは。

    *1:数式載せようかと思ったけど、異常にめんどくさいからやめますごめんなさい。

    【Day-17】DeepLearning系ライブラリ、『Keras』の使い方まとめ(2.x対応版)

    【最終更新 : 2017.12.17】
    ※以前書いた記事がObsoleteになったため、2.xできちんと動くように書き直しました。

    f:id:imslotter:20171217160749p:plain

    データ分析ガチ勉強アドベントカレンダー 17日目。

    16日目に、1からニューラルネットを書きました。 それはそれでデータの流れだとか、活性化関数の働きだとか得るものは多かったのですが、Kerasと言うものを使ってみて、何て素晴らしいんだと感動してしまいました

    今まで苦労して数十行書いていたものが、わずか3行で書ける! 正直、スクラッチで書く意味って、理解にはいいけど研究や分析には必要あんまないんですよね。車輪の再発明になるし。 と言うわけで、使えるものはどんどん使っていこうスタンスで、今日はKerasの紹介です!

    Tutorial+気になった引数を掘り下げて補足のような感じで書いています。 ちなみに、各部のコード以下をつなぎ合わせるとmnistの分類器が動くようになっています。 *1 *2

    Kerasとは?

    Keras is a high-level neural networks library, written in Python and capable of running on top of either TensorFlow or Theano. It was developed with a focus on enabling fast experimentation. Being able to go from idea to result with the least possible delay is key to doing good research. (Documentationより)

    TensorFlowやTheanoのラッパーで、簡単にディープ系の実験環境を整えられるようになっている。 実際に使ってみましたが、かなり直感的にネットワークを構成することが出来る

    使ってみる

    インストール

    前に記事でも少し書いた。

    www.procrasist.com

    pip install tensorflow
    pip install keras

    で入ると思います。その他のライブラリ(opencvとか)は必要に応じて入れてください。

    データを用意

    とりあえずMNIST(デフォルトでmnistデータセットは利用可能) from keras.datasets import hogehogeでhogehogeのデータセットを読み込める。以下は読み込めるデータセット一覧

    データセット 内容 trainデータ数 testデータ数
    mnist 28×28の手書き数字(白黒)10個に分類 60000 10000
    cifar10 32×32のカラー画像10に分類 50000 10000
    cifar100 32×32のカラー画像100に分類 50000 10000
    imdb 25000の映画のレビューのデータセット、ラベルは肯定否定の2種 25000 25000
    reuters 11228個のニュースデータセット46トピックに分類 8972 2246
    fashion_mnist 28×28の白黒画像を10種類のカテゴリに分類 60000 10000

    求めたいものがあれば各々で利用すれば良い。データの渡し方の参考にこれらのデータセットの型を知っておくのがいいかも。

    コード

    from keras.datasets import mnist
    from keras.utils import np_utils
    (X_train, y_train),(X_test,y_test) = mnist.load_data()
    X_train = X_train.reshape(60000,784).astype('float32')
    X_test = X_test.reshape(10000,784).astype('float32')
    #[0,255]の値を[0,1]に正規化
    X_train /= 255.0
    X_test /= 255.0
    # 1 of Kのベクトルに変換
    y_train = np_utils.to_categorical(y_train, 10)
    y_test = np_utils.to_categorical(y_test, 10)
    

    modelを作る

    model = Sequential()で初期化
    model.add()で層を積んでいく

    • 一層目:IN-784次元, OUT-64次元, 活性化関数-Selu
    • 二層目:IN-default(多分64次元) OUT-10次元 活性化関数-SoftMax関数

    入力のイメージは下図。MNISTは28×28の図だが、784次元の1次元ベクトルに変換してから入力している

    • デフォである活性化関数
    • より高度な活性化関数 : keras.layers.advanced_activations モジュール内

    関数の違いは下図 f:id:imslotter:20170107153956p:plain

    コード

    from keras.models import Sequential
    from keras.layers import Dense
    model = Sequential()
    model.add(Dense( activation="selu", units=64,input_dim=784))
    model.add(Dense(activation="softmax", units=10))
    

    学習プロセス

    model.compile

    • 損失関数 : どのくらいの誤差があるかを定量
    • 最適化関数 : 誤差からのパラメータ更新の仕方

    loss : 損失関数

    デフォで用意されている関数

    関数名 意味
    mean_squared_error, mse 二乗誤差
    mean_absolute_error, mae 絶対誤差
    mean_absolute_percentage_error, make 正解とのズレ(絶対値)の割合
    mean_squared_logarithmic_error, msle 正解とのズレ(二乗誤差)の割合
    squared_hinge マイナスのところは0, プラスは二乗誤差
    hinge マイナスのところは0, プラスは絶対値
    binary_crossentropy loglossとしても知られている
    categorical_crossentropy マルチクラスloglossとして知られている(ラベルがバイナリ配列であることが必要)
    sparse_categorical_crossentropy スパースラベルを取る
    kullback_leibler_divergence, kld 確率を分布とみなした時の差
    poisson 予測-正解*log(予測)の平均
    cosine_proximity 予測と正解間のコサイン近似の負の平均

    詳しい実装例はこちら GitHub

    optimizer : 最適化

    デフォで用意されている関数

    • SGD
    • RMSprop(デフォルトパラメータ推奨)
    • Adagrad(デフォルトパラメータ推奨)
    • Adadelta(デフォルトパラメータ推奨)
    • Adam(デフォルトは提案論文通り)
    • Adamax Adamの拡張(無限ノルム)
    • Nadam Adam+RMSprop with momentumらしい
    • TFOptimizer よくわかんない
    • 引数の詳細な設計も可能(こちらを参考 web)
    • 新しい最適化関数を作りたいなら GitHubを参考に

    指標の参考には以前書いたこちらの記事もどうぞ

    コード

    from keras.optimizers import SGD
    # optimizer = "sgd"でも、簡単に最適化関数の設定ができる
    #(その場合パラメータはデフォルト値)
    model.compile(loss="categorical_crossentropy", 
                  optimizer=SGD(lr=0.01, momentum=0.9, nesterov=True),
                  metrics=["accuracy"])
    

    モデルの学習

    学習

    model.fit()でモデルを学習。引数は以下

    • X_train : データ(MNISTだと、どこに何個データが入っているか)
    • y_train : ラベル(1-10)
    • batch_size : ミニバッチにデータを分けるこの場合60000個ある学習データを、32個ずつに分けて、重みを更新する
    • epochs : 学習の繰り返し。60000個を何回学習するか
    • validation_split : 学習データの何パーセントをvalidation用データにするか(各エポックのテスト。validationの説明はこの記事を参考にどうぞ : 【Day-10】Cross Validationとパラメータサーチでモデルの調整 )

    また、学習の様子も保存されている model.fit内の変数を見てみるとこんな感じ

    print(history.__dict__)
    >>> {'model': <keras.models.Sequential object at 0x110c78510>, 'params': {'verbose': 1, 'nb_epoch': 3, 'batch_size': 32, 'metrics': ['loss', 'acc', 'val_loss', 'val_acc'], 'nb_sample': 48000, 'do_validation': True}, 'epoch': [0, 1, 2], 'history': {'acc': [0.89862500000000001, 0.94735416666666672, 0.96150000000000002], 'loss': [0.35020409799863894, 0.18034435256198048, 0.132605804043822], 'val_acc': [0.94125000000000003, 0.95816666666666672, 0.96499999999999997], 'val_loss': [0.20651897403846184, 0.14732480475058157, 0.1263637450747192]}}
    

    コールバック関数

    keras.callbacks内の関数を使ってモデルの保存ができる。

    • check = ModelCheckpoint("modelname.hdf5")で、学習済みモデルの重みを保存する
      • その場合、model.fit(hogehoge,callbacks=[check])としてやる
    • tensorboardなどとも連携可能
    • 自分でcallbackしたい関数を作ることも可能
    • 詳しくはドキュメント参照

    コード

    from keras.callbacks import ModelCheckpoint
    check = ModelCheckpoint("model.hdf5")
    history = model.fit(X_train, y_train, epochs=20, 
                        validation_split=0.2, batch_size=32,
                        callbacks=[check])
    

    モデルの評価

    model.evaluate()という関数で、テストデータを用いたモデルの評価が可能。lossとaccuracyを見ている

    コード

    loss, accuracy = model.evaluate(X_test, y_test)
    print("\nloss:{} accuracy:{}".format(loss, accuracy))
    

    モデルの可視化

    可視化(tensorflow)

    モデルの可視化

    ニューラルネットのネットワークを可視化するとき
    keras.utils.visualize_util
    手順(Mac, Homebrew有り)

    • $brew install graphviz
    • $pip install pydot

    とすれば使える

    • to_file : 出力の画像ファイルの名前
    • show_shapes : グラフ中に出力形状を書くか(デフォはFalse)
    • show_layer_names : レイヤー名を書くか(デフォはTrue)

    コード

    from keras.utils.visualize_util import plot
    plot(model, to_file="model.png", show_shapes=True, show_layer_names=True)
    

    こんなのが出てくる
    f:id:imslotter:20170107113053p:plain

    学習の様子をplot

    matplotlibで可視化

    コード

    import matplotlib.pyplot as plt
    
    def plot_history(history):
        # 精度の履歴をプロット
        plt.plot(history.history['acc'],"o-",label="accuracy")
        plt.plot(history.history['val_acc'],"o-",label="val_acc")
        plt.title('model accuracy')
        plt.xlabel('epoch')
        plt.ylabel('accuracy')
        plt.legend(loc="lower right")
        plt.show()
    
        # 損失の履歴をプロット
        plt.plot(history.history['loss'],"o-",label="loss",)
        plt.plot(history.history['val_loss'],"o-",label="val_loss")
        plt.title('model loss')
        plt.xlabel('epoch')
        plt.ylabel('loss')
        plt.legend(loc='lower right')
        plt.show()
    # modelに学習させた時の変化の様子をplot
    plot_history(history)
    

    こんな感じに

    f:id:imslotter:20170107151009p:plain

    正解率98%とかって半端ないですね。人間と同じかそれ以上の識別性能かも!

    Keras 1.xからKeras 2.xの変更点

    詳しくはこちらの公式記事を参考

    上記チュートリアルで変更した点をメモとして記しておく。

    Layerを作る際に、活性化関数も一緒に入れるようになった

    from keras.layers import Dense, Activation
    model = Sequential()
    model.add(Dense(output_dim=64, input_dim=784))
    model.add(Activation("relu"))
    

    出力

    UserWarning: Update your `Dense` call to the Keras 2 API: `Dense(input_dim=784, units=64)

    変更後

    model = Sequential()
    model.add(Dense( activation="relu", units=64,input_dim=784))

    nb_epoch -> epochs

    from keras.callbacks import ModelCheckpoint
    check = ModelCheckpoint("model.hdf5")
    history = model.fit(X_train, y_train, nb_epoch=20, 
                        validation_split=0.2, batch_size=32,
                        callbacks=[check])
    

    出力

    UserWarning: The `nb_epoch` argument in `fit` has been renamed `epochs`.
      warnings.warn('The `nb_epoch` argument in `fit` '

    変更後

    from keras.callbacks import ModelCheckpoint
    check = ModelCheckpoint("model.hdf5")
    history = model.fit(X_train, y_train, epochs=20, 
                        validation_split=0.2, batch_size=32,
                        callbacks=[check])
    

    show_accuracyの廃止 : defalutでaccuracyが出力されるようになった。

    model.evaluate(X_test, y_test, show_accuracy=True)
    

    出力

    TypeError: evaluate() got an unexpected keyword argument 'show_accuracy'

    変更後

    model.evaluate(X_test, y_test)
    

    まとめ

    今日のポイントは
    ・Kerasを使うと、簡単にディープラーニングができる
    ・可視化も簡単にできるので、実験にはもってこい
    ・ただ、便利・簡単すぎるので、中身の仕組み知りたい人は一から書いてみるのもいいかも!

    実際にモデルを動かしてる部分って、たかだか十数行とかですものね。すごい。 いろいろな実験がしやすそうだし、これからどんどん使っていきたい。

    なお、Kerasを書籍で勉強したい方はこちらがオススメです。特にRNNに関しては丁寧に書かれているので!

    詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

    詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~

    明日はKerasを用いたRNN実装にチャレンジしてみる(予定)です!お楽しみに!

    *1:日本語があると怒られる方は、ファイル先頭に#coding:utf-8をつけておきましょう。

    *2:順を追ってimportしているので、そのまま使うとものすごく見にくいかもです。コピペした後は適当に体裁を整えてください

    【Day-16】ニューラルネットを0から作り、仕組みを基礎から理解する

    f:id:imslotter:20171216170643p:plain

    データ分析ガチ勉強アドベントカレンダー 16日目。

    今日からは少しディープラーニングの勉強。

    ここ数年間、深層学習用ライブラリも猛烈に整備され、誰でも簡単にディープラーニングを使えるようになりました。

    その一方で、整備されすぎて、魔法の箱だという認識も多いですよね。 けれど、深層学習と言えど、しているのはほとんど線形代数微積分を組み合わせた数値計算です。

    だったら自分で作れるのでは? というわけで、仕組みを理解するために、0からスクラッチで作ることにしました。

    尚、勉強にはプロフェッショナルシリーズの深層学習を利用しています。

    爆速で技術が進む深層学習界隈では少々obsoleteかもしれませんが*1、きちんと基礎の基礎を知るにはいい本だと思います。詳しい計算方法を学びたい人は、どうぞ。(線形代数偏微分の知識が必要です。)

    深層学習 (機械学習プロフェッショナルシリーズ)

    深層学習 (機械学習プロフェッショナルシリーズ)

    作るニューラルネットワーク

    下図のような基本的な3層ニューラルネットワーク

    • 全結合
    • 活性化関数 : 隠れ層 : Sigmoid関数, 出力層 : ソフトマックス関数
    • ミニバッチ処理を採用

    の構成で作ってみる

    ニューラルネットの設計

    ニューラルネットには下記の要素が必要。

    • ネットワーク構成 : 何段構成にするとだとか、どういう計算を施すか解か、どこの層をつなぐとかを考えなければならない。今回は3段の単純なネットワークだけど、実際は用途によって様々な構成が出来る。Network Zooなどが、ネットワーク構成の参考になる。けど、今でも試行錯誤

    • 活性化関数 : ニューラルネット非線形性を生み出す重要な要素。

    • 重みのアップデート方法 : ニューラルネット学習を進める上で重要な要素。

    活性化関数も、試行錯誤で作られている事が多い。 近年よく使われるのはReLU等。活性化関数に関しては、以前記事にしたことがあるのでそちらを参考に。

    www.procrasist.com

    順伝搬計算

    入力(ベクトル)から、線形変換と活性化関数による非線形計算をしながら出力層の値を取り出す。 最終的に各値のスコアが出てくるので、回帰ならそのまま用いればいいし、分類ならスコアが最も高いものを選ぶ。

    逆伝搬計算

    出力値から、誤差を出すことが出来る。その誤差を元にして、各層の重みを調整していく。 このためには微分の知識が必要になってくる。出力層から入力層計算を進めるため、誤差逆伝搬と呼ばれる。 誤差にはクロスエントロピー関数を用いる事が多い

    重みの更新

    逆誤差伝搬を行うことで、誤差からの各重みが出せる。最もシンプルなのは、単に誤差を一回の微分値で調整する方法(勾配降下法)。確率的勾配法(SGD)などもこれを元にしている。 この更新方法も様々に考えられていて、もう少し高度なモメンタムや、学習率を自動的に更新してくれるAdam, AdaGradなどもある。

    f:id:imslotter:20171216162434p:plain
    (重みの更新式{\epsilon}は学習率)

    実装

    以上を踏まえて、実装をしてみた。いつものように

    ネットワークの設計

    必要なネットワークの構成と、活性化関数、順伝搬計算、逆伝搬計算を、実装

    コードをみる
    class NN:
        def __init__(self, num_input, num_hidden, num_output, learning_rate):
            self.num_input = num_input
            self.num_hidden = num_hidden
            self.num_output = num_output
            self.learning_rate = learning_rate
    
            self.w_input2hidden = np.random.random((self.num_hidden, self.num_input))
            self.w_hidden2output = np.random.random((self.num_output, self.num_hidden))
            self.b_input2hidden = np.ones((self.num_hidden))
            self.b_hidden2output = np.ones((self.num_output))
    
        ##活性化関数(シグモイド関数)
        def activate_func(self, x):
            return 1/(1+np.exp(-x))
        ##活性化関数の微分
        def dactivate_func(self,x):
            return self.activate_func(x)*(1-self.activate_func(x))
        ##ソフトマックス関数
        def softmax_func(self,x):
            C = x.max()
            f = np.exp(x-C)/np.exp(x-C).sum()
            return f
        ##順伝播計算
        def forward_propagation(self, x):
            u_hidden = np.dot(self.w_input2hidden, x) + self.b_input2hidden
            z_hidden = self.activate_func(u_hidden)
            u_output = np.dot(self.w_hidden2output, z_hidden) + self.b_hidden2output
            z_output = self.softmax_func(u_output)
            return u_hidden, u_output, z_hidden, z_output
        ##逆伝播でdeltaを求める
        def backward_propagation(self,t,u_hidden,z_output):
            t_vec = np.zeros(len(z_output))
            t_vec[t] = 1
            delta_output = z_output - t_vec
            delta_hidden = np.dot(delta_output, self.w_hidden2output * self.dactivate_func(u_hidden))
            return delta_hidden, delta_output
        ##wに関するgradient
        def calc_gradient(self,delta,z):
            dW = np.zeros((len(delta), len(z)))
            for i in range(len(delta)):
                for j in range(len(z)):
                    dW[i][j] = delta[i] * z[j]
            return dW
        # update(SGD)
        def update_weight(self,w0,gradE):
            return w0 - self.learning_rate*gradE
    

    学習

    ランダムにデータを選ぶことで、確率的勾配法にしてみた。

    コードをみる
    def train(nn, iteration,savefig=False):
        epoch = 0
        for epoch in range(iteration+1):
            grad_i2h = 0
            grad_h2o = 0
            gradbias_i2h = 0
            gradbias_h2o = 0
            rand = randint(0,len(data),100)
            for r in rand:
                u_hidden, u_output, z_hidden, z_output = nn.forward_propagation(data[r])
                delta_hidden, delta_output = nn.backward_propagation(target[r], u_hidden, z_output)
                grad_i2h += nn.calc_gradient(delta_hidden, data[r])
                grad_h2o += nn.calc_gradient(delta_output, z_hidden)
                gradbias_i2h += delta_hidden
                gradbias_h2o += delta_output
            nn.w_input2hidden = nn.update_weight(nn.w_input2hidden, grad_i2h / len(rand))
            nn.w_hidden2output = nn.update_weight(nn.w_hidden2output, grad_h2o / len(rand))
            nn.b_input2hidden = nn.update_weight(nn.b_input2hidden, gradbias_i2h / len(rand))
            nn.b_hidden2output = nn.update_weight(nn.b_hidden2output, gradbias_h2o / len(rand))
    

    実行

    sklearnでデータセットを用意し、実際に分類できるかを分離平面で表してみている。 自信度によって色の濃さを調整しているので、学習が進んでいくと、きちんと分類されて、濃くなっていくはず。。。!

    コードをみる
    #データの用意
    num_cls = 5
    data,target = make_blobs(n_samples=1000, n_features=2, centers=num_cls)
    # Neural Netの用意
    nn = NN(num_input=2,num_hidden=20,num_output=num_cls,learning_rate=0.2)
    # 学習
    train(nn, iteration=500, savefig=True)
    
    plt.figure(figsize=(5,5))
    #色の用意
    base_color = ["red","blue","green","yellow","cyan",
                  "pink","brown","gray","purple","orange"]
    colors = [base_color[label] for label in target]
    # 教師データのプロット
    plt.scatter(data[:,0],data[:,1],color=colors,alpha=0.5)
    print("plotting...")
    xx = np.linspace(-15,15,100)
    yy = np.linspace(-15,15,100)
    for xi in xx:
        for yi in yy:
            _,_,_,z_output = nn.forward_propagation((xi, yi))
            cls = np.argmax(z_output) #softmaxのスコアの最大のインデックス
            score = np.max(z_output)
            plt.plot(xi, yi, base_color[cls],marker="x",alpha=s)
        print(".",end="")
    plt.xlim=(-15,15)
    plt.ylim=(-15,15)
    plt.show()
    print("finish plotting")
    

    実行結果

    以下が、300イテレーション回したときの分離平面の様子である。

    また、学習が進んでいく様子もアニメーションにしてみた。(画像が粗くてすみません)

    うまく学習が進んでいる様子が見て取れます。

    まとめ

    今回は実際に作ることで、ニューラルネットの仕組みを理解してみた。 実際に自分で作るとどこで非線形な要素が生まれるのかとか、誤差をどう学習に反映させているのかというのが分かるようになるので、勉強になります。魔法の箱なんかじゃなく、ちゃんとした学習器として理解できたように思います。0空作ればここから柔軟にいろいろと発展させていくことも可能ですしね。

    今回は式からプログラムを書き起こしたが、最近は↓のようなスクラッチで書く人のためのとても良い本も売っているので、そちらも参考にしてみてはいかがでしょうか。

    また、この辺の基礎を知っている人は、正直ライブラリを使ったほうが早いです。 明日からは、深層学習用ライブラリをいくつか触っていきたいと思います。(KerasかPytorch) ではでは!

    *1:言うても2-3年前までは最先端

    【Day-15】ベイズ的最適化で最強のゴールデンクロスを見つける

    f:id:imslotter:20171215132801p:plain

    データ分析ガチ勉強アドベントカレンダー 15日目。

    最強のゴールデンクロス、それは最も儲かるように移動平均線を引いたときの交点 である。

    一説には、テクニカル分析は、チャートにすべての情報が詰まっているという前提があるそうですね。 ということは、データから最適な移動平均を求められるのでは?

    というわけで、今回はベイズ最適化の勉強もかねて、最適な移動平均のwindow幅をデータから出してみようと思います。

    ※始めに実験を載せてます。理論が気になる方は理論のところから見ると良いでしょう。 [:contents]

    ゴールデンクロス/デッドクロスとは?

    昨日の記事をちょっとだけ復習。

    www.procrasist.com

  • ゴールデンクロス : 株価が下落した後に短期移動平均線が、長期移動平均線下から上に抜ける現象。買いのサイン
  • デッドクロス : 株価が上昇した後に短期移動平均線が、長期移動平均線上から下に抜ける現象。売りのサイン
  • (わかる株式用語より)

    長期/短期2本の移動平均線が必要なわけです。

    【実験】

    実験設定

    二つの移動平均線を使って、売り買いを判断するシステムがあると仮定する。そいつの動きはこう

    • 2015年1月4日、資産100株&100万円からスタート
    • データとして使うのは3年分の日経平均
    • ゴールデンクロス買い(全資産)デッドクロス売り(全株)
    • 2017/12/13時点で株を全部売って、総資産を計算
    • 評価値は、何%資産が増えたか*1

    つまり、シグナルにあわせて売り注文、買い注文をしていき、資産の増分が評価になる。 昨日作ったアルゴリズムで、ゴールデンクロスデッドクロスを発見できるようになっているので、そこに売り買いを入れるだけ。いいタイミングで売り買いをし、波を捉えられるかというゲーム。

    最適化

    最適化については最後に理論的な説明を簡単につけているので、そこを参照。ベイズ的最適化を用いる。下記がその時に扱うものの設定

    要素 選択
    Acquisition function EI (Acqusition functionの参考)
    kernel Mattern52
    イテレーション回数 50

    尚、実装はgithubにまとめている。

    結果は・・・?

    さあ、ベイズ最適化で計算してみた結果どうなったか!! まずは図を御覧ください。

    最適化してみた結果の、平均、標準偏差、acquisition funciton *2

    試行を繰り返せば、一番右の図の濃い黄色の部分が、選ぶべき最適な点となっているはずです。

    この結果、

    window1=91, window2 = 208で移動平均を取ると、最も儲かる

    という結果を得ることができました。ちなみにこのとき、50%資産増でした。3ヶ月と7ヶ月で移動平均ですね!!

    ちなみに、その場合に移動平均を引いたらこうなる。

    あんまり売買しないほうが儲かるということかな?ちなみにまだ売りのタイミングは来ていないようです。 ここでのクロスポイントで、売るといいのかな?データを分析したらそういってる!!

    【実装】

    GPyOpt/GPy

    今回は、実装にGPyOptを利用することにする。GPyOptは、ガウス過程計算ライブラリGPyを利用してベイズ的最適化を行うライブラリ。なので、どちらもインストールしておく必要がある。

    pip install -U scipy #scipyも最新版にしておく
    pip install GPy
    pip install gpyopt

    ※実装は長くなるので、記事では割愛しますgithubを参考にしてください。評価関数を作って、GPyOptに投げ入れているだけです!(day15)

    github.com

    【理論】ベイズ的最適化とは?

    荒野の中から宝を探す

    こちらの記事が分かりやすかったです。ガウス過程に従うサンプルとしてモデル化して最適化を行うような手法

    難しいことは下記の参考に譲るとしましょう。どれも分かりやすくていい資料です!

    ざっくりとポイントだけ。

  • f(x)の最大値を求めたい
  • fの形は分かっていない
  • けれど、xを与えて、f(x)を計算するまでに時間がかかる
  • f(x)を何か別の関数で近似する方法がほしい
  • 普通の式ならいいのですが、たとえば、xがパラメータ、f(x)が何かしらのシミュレータの場合、 まずf(x)の形を定式化することが難しいです。コレをブラックボックスといいます。さらに、一回一回の評価にめちゃくちゃ時間がかかる**とすると、xをちょっとずつ変えて全探索など、厳しいですね。

    たとえて言うなら、荒野の中から、当てもなく土を掘り返して中に眠る宝を探すイメージです。途方もないですね。 その絶望的な状況でも、賢く範囲を絞っていこうというのが、ベイズ最適化です。

    関数近似

    ベイズ最適化を用いると、平均と標準偏差の確率分布が出てくる。点が与えられた場合はその値に収束するので自信を持ってそうだといえるが、それ以外の値はあてずっぽうなので、点以外の場所は自身がなくなる(≒分散が大きくなる)。つまり、グラフは下記のようになる。



    A Tutorial on Bayesian Optimization ...より。新しい点が与えられると、自信のない領域(青くマスクされた領域)が減っていることが確認できる。

    最適化のための指針、Acquisition function

    また、上の図のAcquisition Functionというのが、最適化のために次の点を探す指針となる。 acquisition関数にはいろんな物を設定できるが、EIやUCBなどが代表的。UCBの式は、平均と標準偏差の確率分布で、下記のようにあらわされる

    \displaystyle{UCB=\mu(x)+\sqrt{\beta}\sigma(x)}

    この関数の最大値を次の点に選ぶ。 第一項目の{\mu}が近似関数なのでその最大値を選べばよいが、自信のないところも探したいそれが第二項である*3。ここから新しい点を選んで、また関数近似してを繰り返すのが、ベイズ的最適化です。

    注意

    とても強力な手法だが、下記あたりはアルゴリズム上の欠点として注意したいところ。

    • 関数近似逆行列計算が入るので(データ数)3オーダーで計算が重くなっていく
    • どこからも離れた点では分散が大きくなるので、次元が大きくなると**選ばれる点が端っこに寄る

    ※今ではこれを解決するようないろんな工夫もなされているっぽい。

    まとめ

    実験的に最適化を行ってみました。ここ三年程好景気なので、波をつかめばかなり儲かるっぽいですね。あとはチャートと自分の作ったアルゴリズムを信じるのみ!!

    また、ベイズ的最適化は、探索と活用のバランス次第では、局所解に陥る可能性もあるので、試行によっては解が変わることもあります。安定させるためにまだ色々と改良できると思っていますので、試行錯誤を繰り返していきます。もうちょっとテストもしないといけなさそうだなぁ。

    とりあえず時系列データはこのあたりで一旦置いておいて、明日は深層学習について触れたいと思っています。それではまた明日!

    *1:手数料は考えない

    *2:実装上最小値にしたかったので、評価値に-1をかけている

    *3:いわゆる探索と活用のトレードオフ

    【Day-14】株価や仮想通貨で使える、5つのテクニカル分析を解説&Pythonで実装してみた

    f:id:imslotter:20171214165640p:plain

    データ分析ガチ勉強アドベントカレンダー 14日目。

    時系列データでまず思いつくのは、株価のチャートですよね。 また、最近はやっている仮想通貨。私も最近coincheckに入金しました。

    やっぱ、実際にお金が絡むとちゃんと勉強しようって言う気になる!笑

    せっかくチャートを見るわけだし、その見方について勉強しておこうと思いました。 そしてせっかくなので、自分で実装してどういう仕組みなのかまで知っておこうと思いました。 理系だからね、分からないものを使うのは嫌だからね。

    というわけで、Python(主にPandasとMatplotlibを用いながら)でテクニカル指標についてやっていきます。扱うデータは三年分の日経平均株価

    指標について知りたい人も、自分で実装してみたいという人もどうぞ。

    テクニカル分析ファンダメンタル分析

    大きく相場にはテクニカル分析ファンダメンタル分析の2つの分析手段がある。

    • ファンダメンタル分析 : 情報を集め、価格がどうなるかを判断する。株で言えば会社の情報とか、仮想通貨で言えば技術的な要素とか方針とかそういうのかな?
    • テクニカル分析 : チャート情報から、価格がどうなるかを判断する。チャートだけ見れば良いが、大変動などの予測は難しい。

    一般的に、テクニカルのほうがとっつきやすいけど、判断を誤ることも多いと言われている。

    トレーダーが良く使うテクニカル指標より、主要なテクニカル分析の指標を並べてみた。

    順位 テクニカル指標 割合
    1位 移動平均線 36.6%
    2位 ボリンジャーバンド 17.1%
    3位 MACD 13.2%
    4位 RSI 8.4%
    5位 一目均衡表 8.4%
    6位 ストキャスティクス 7.1%
    7位 フィボナッチ 2.9%
    8位 その他 6.3%

    これの、1位から5位までの指標について簡単に説明&Python実装する。

    実装において

    実装を見る方は、まず下記でデータを読み込んで置いてください。また、自分の好きなデータでもOKです。

    いつものごとく、githubのノートブックに上げています。(day14)

    github.com

    コードをみる
    # データの読み込み
    %matplotlib inline
    import pandas as pd
    import numpy as np
    import matplotlib.pyplot as plt
    from datetime import datetime
    def str2time(string):
        return datetime.strptime(string, '%Y/%m/%d')
    df_orig = pd.read_csv("stock_nikkei_3yrs.csv")
    df_orig.index = df_orig.date.apply(str2time)
    df_orig = df_orig.sort_index(ascending=True)
    

    ローソク足

    説明

    チャートの代表。一つ一つのローソクを見ると、ひと目でその期間の動きが分かるようになっている。

  • 始値・高値・安値・終値をまとめて表現
  • 終わりのほうが始まりより高いと青色 。逆だと赤色
  • 上下に伸びた線はその期間の最高値と最低値。ヒゲという
  • https://www.fxbroadnet.com/techs/rousokuahi.gif
    テクニカルチャート講座:ローソク足より引用

    実装

    コードをみる
    import matplotlib.finance as mpf
    from matplotlib.dates import date2num
    df = df_orig[df_orig.index >= datetime(2017,7,1)]
    df_ohlc = pd.DataFrame(df[["open","high","low","close"]])
    xdate = [x.date() for x in df_ohlc.index] #Timestamp -> datetime
    ohlc = np.vstack((date2num(xdate), df_ohlc.values.T)).T #datetime -> float
    
    fig = plt.figure(figsize=(15,5))
    ax = plt.subplot()
    mpf.candlestick_ohlc(ax, ohlc, width=0.7, colorup='b', colordown='r')
    ax.grid() #グリッド表示
    ax.set_xlim(xdate[-60], xdate[-1]) #x軸の範囲
    fig.autofmt_xdate() #x軸のオートフォーマット
    

    1位 : 移動平均線

    説明

    そのままのチャートはノイズ等も多く、価格が乱高下するので、その日からある期間分の平均を取るという操作を良くする。英語で言うとMoving Averageで、25MAだと25日分の移動平均線ということになる。平均を取ると、ノイズが無くなり、そのチャートのトレンドが明確に現れる。

    また、2つの移動平均線を組み合わせることも多く、各線の交差するところをゴールデンクロス/デッドクロスと呼ぶ。

  • ゴールデンクロス : 株価が下落した後に短期移動平均線が、長期移動平均線下から上に抜ける現象。買いのサイン
  • デッドクロス : 株価が上昇した後に短期移動平均線が、長期移動平均線上から下に抜ける現象。売りのサイン
  • (わかる株式用語より)
    • 日足の場合、短期が25日、長期が75日
    • 週足の場合、短期が13週、長期が26週

    が基本らしい。

    実装

    コードをみる
    df1 = df_orig.sort_index()
    df1["ma25"] = df1.close.rolling(window=25).mean()
    df1["ma75"] = df1.close.rolling(window=75).mean()
    df1["diff"] = df1.ma25-df1.ma75
    # df1 = df.date.apply(datetime.timestamp())
    df1["unixtime"] = [datetime.timestamp(t) for t in df_orig.index]
    plt.figure(figsize=(15,5))
    
    
    xdate = [x.date() for x in df1.index]
    # line and Moving Average
    plt.plot(xdate, df1.close,label="original")
    plt.plot(xdate, df1.ma75,label="75days")
    plt.plot(xdate, df1.ma25,label="25days")
    plt.xlim(xdate[0],xdate[-1])
    plt.grid()
    
    #Cross
    for i in range(1, len(df1)):
        if df1.iloc[i-1]["diff"] < 0 and df1.iloc[i]["diff"] > 0:
            print("{}:GOLDEN CROSS".format(df1.iloc[i]["date"]))
            plt.scatter(df1.iloc[i]["date"],df1.iloc[i]["ma25"],marker="o",s=100,color="b")
            plt.scatter(df1.iloc[i]["date"],df1.iloc[i]["close"],marker="o",s=50,color="b",alpha=0.5)
    
        if df1.iloc[i-1]["diff"] > 0 and df1.iloc[i]["diff"] < 0:
            print("{}:DEAD CROSS".format(df1.iloc[i]["date"]))
            plt.scatter(df1.iloc[i]["date"],df1.iloc[i]["ma25"],marker="o",s=100,color="r")
            plt.scatter(df1.iloc[i]["date"],df1.iloc[i]["close"],marker="o",s=50,color="r",alpha=0.5)
    plt.legend()
    

    2位 : ボリンジャーバンド

    説明

    移動平均線移動標準偏差を用いて分析を行う。

    {1\sigma, 2\sigma, 3\sigma}線を移動平均に引く

    • {1\sigma} : 68%の確率で株価がそこに収まる
    • {2\sigma} : 95%の確率で株価がそこに収まる
    • {3\sigma} : 99%の確率で株価がそこに収まる

    均衡状態の相場での正規分布を仮定しているので、 ここから外れた = 釣り合いが破れているということ。

    不安定だが、売買のチャンスともいえる。

    株の売買サインを簡単に見極めるボリンジャーバンドの使い方によると

  • バンドが狭いときは様子見(スクイーズ状態)
  • バンドが上下に開いているときはトレンドに入り気味。上昇のときはチャンス
  • +2σのラインを価格がぶち破ればめちゃくちゃ強いトレンド
  • 実装

    コードをみる
    def Bollinger(df,window=25):
        df1 = df.copy()
        df1["ma"] = df1.close.rolling(window=window).mean()
        df1["sigma"] =  df1.close.rolling(window=window).std()
        df1["ma+2sigma"] = df1.ma + 2*df1.sigma
        df1["ma-2sigma"] = df1.ma - 2*df1.sigma
        df1["diffplus"] = df1.close - df1["ma+2sigma"] 
        df1["diffminus"] = df1["ma-2sigma"] - df1.close
        s_up = df1[df1["diffplus"]>0]["close"]
        s_down = df1[df1["diffminus"]>0]["close"]
    
        plt.figure(figsize=(15,5))
        plt.grid()
        plt.scatter(s_up.index, s_up.values,marker="x",s=100,color="blue")
        plt.scatter(s_down.index, s_down.values,marker="x",s=100,color="red")
        xdate = [x.date() for x in df1.index]
        plt.plot(xdate, df1.close.values,label="close",color="b",alpha=0.8)
        plt.plot(xdate,df1.ma.values,label="{}ma".format(window))
        plt.fill_between(xdate, df1.ma-df1.sigma,df1.ma+df1.sigma,color="red", alpha=0.7, label="$1\sigma$")
        plt.fill_between(xdate, df1.ma-2*df1.sigma,df1.ma+2*df1.sigma,color="red", alpha=0.3, label="$2\sigma$")
    #     plt.fill_between(xdate, df1.ma-3*df1.sigma,df1.ma+3*df1.sigma,color="red", alpha=0.3, label="$3\sigma$")    
        plt.legend()
        plt.xlim(xdate[0],xdate[-1])
    Bollinger(df_orig,window=25)
    

    3位 : MACD(Moving Average Convergence/Divergence)

    説明

    移動平均線を少し改良している。時系列は直近のデータのほうが昔のデータよりも価値が高い傾向にあるので、そういう平均のとり方をして、移動平均線を引いている(指数平滑移動平均線(EMA)という)。下記の2本の移動平均線を見ることで、そのトレンドを把握する。

    • MACD : 12日EMA - 26日EMA
    • シグナル : MACDの9日EMA
  • シグナルをMACD下から上へクロス。買いのサイン
  • シグナルをMACD上から下へクロス。売りのサイン
  • 実装

    コードをみる
    df_orig.head()
    df_orig["MACD"] = df_orig.close.ewm(span=12, min_periods=1).mean() - df_orig.close.ewm(span=26, min_periods=1).mean()
    df_orig["signal"] = df_orig.MACD.ewm(span=9, min_periods=1).mean()
    df_orig["macd_diff"] = df_orig["MACD"]-df_orig["signal"]
    xdate = [x.date() for x in df_orig.index]
    
    plt.figure(figsize=(15,10))
    plt.subplot(211)
    plt.plot(xdate, df_orig.close,label="original")
    plt.xlim(xdate[0],xdate[-1])
    plt.legend()
    plt.grid()
    plt.subplot(212)
    plt.title("MACD")
    plt.plot(xdate, df_orig.MACD,label="MACD")
    plt.plot(xdate, df_orig.signal,label="signal")
    plt.xlim(xdate[0],xdate[-1])
    plt.legend()
    plt.grid(True)
    #Cross
    for i in range(1, len(df_orig)):
        if df_orig.iloc[i-1]["macd_diff"] < 0 and df_orig.iloc[i]["macd_diff"] > 0:
            print("{}:GOLDEN CROSS".format(df_orig.iloc[i]["date"]))
            plt.scatter(xdate[i],df_orig.iloc[i]["MACD"],marker="o",s=100,color="b")
    
        if df_orig.iloc[i-1]["macd_diff"] > 0 and df_orig.iloc[i]["macd_diff"] < 0:
            print("{}:DEAD CROSS".format(df_orig.iloc[i]["date"]))
            plt.scatter(xdate[i],df_orig.iloc[i]["MACD"],marker="o",s=100,color="r")
    

    4位 : RSI (Relative Strength Index)

    買われすぎ売られすぎな相場では、投資家が不安になってくる。その心理状況を定量化しようとした指標。 値上がりと値下がりのバランスを見ている。

  • RSIの数値が70%以上になると買われ過ぎゾーン
  • RSIの数値が30%以下になると売られ過ぎゾーン
  • それぞれのゾーンに入った後に反転した動きになったところで買われ過ぎゾーンの時は“売り”を、売られ過ぎゾーンの時は“買い”を考える。
  • カブドットコム証券より

    実装

    コードをみる
    def plot_RSI(df,window):
        diff = df_orig.close.diff(periods=1).values
        xdate = [x.date() for x in df_orig.index]
        RSI = []
        for i in range(window+1,len(xdate)):
            neg = 0
            pos = 0
            for value in diff[i-window:i+1]:
                if value > 0:
                    pos += value
                if value < 0:
                    neg += value
            pos_ave = pos/window
            neg_ave = np.abs(neg/window)
            rsi = pos_ave/(pos_ave+neg_ave)*100
            RSI.append(rsi)
        #RSIのグラフ描画
        plt.plot(xdate[window+1:], RSI,label="RSI {}".format(window),lw=2.5,alpha=0.6)
        plt.xlim(xdate[window+1],xdate[-1])
        plt.ylim(0,100)
        plt.legend()
    xdate = [x.date() for x in df_orig.index]
    plt.figure(figsize=(15,10))
    plt.subplot(211)
    plt.plot(xdate, df_orig.close,label="original")
    plt.xlim(xdate[0],xdate[-1])
    plt.legend()
    plt.grid()
    plt.subplot(212)
    plt.grid()
    plt.title("RSI")
    plot_RSI(df_orig,window=14)
    plot_RSI(df_orig,window=28)
    plot_RSI(df_orig,window=105)
    plt.fill_between(xdate,np.ones(len(xdate))*30,color="blue",alpha=0.2)
    plt.fill_between(xdate,np.ones(len(xdate))*70,np.ones(len(xdate))*100,color="red",alpha=0.2)
    plt.plot(xdate,np.ones(len(xdate))*30,color="blue",linestyle="dotted")
    plt.plot(xdate,np.ones(len(xdate))*70,color="red",linestyle="dotted")
    plt.show()
    

    5位 : 一目均衡表

    説明

    これもローソク足と並ぶ、国産のテクニカル分析方法

    性質
    転換線(Conversion Line) (9日間の高値+9日間の安値)÷2
    基準線(Base Line) (26日間の高値+26日間の安値)÷2
    先行スパン1(Leading Span 1) (転換線+基準線)÷2を26日先にプロット
    先行スパン2(Leading Span 2 (52日間の高値+52日間の安値)÷2を26日先にプロット
    遅行スパン(Lagging Span) 当日の終値を26日前にプロット
  • 転換線と基準線の交差ポイントは、トレンドの傾向。転換線が下から上好転上から下逆転
  • 先行スパン1と先行スパン2の間をという。株価が雲よりも上にいるときは上昇基調、株価が雲より下にいるときは下落基調
  • 実装

    コードをみる
    df_orig = pd.read_csv("stock_nikkei_3yrs.csv")
    df_orig.index = df_orig.date.apply(str2time)
    df_orig = df_orig.sort_index(ascending=True)
    
    max_9 = df_orig.high.rolling(window=9).max()
    min_9 = df_orig.high.rolling(window=9).min()
    df_orig["tenkan"] = (max_9+min_9)/2
    df_orig["base"] = (df_orig.high.rolling(window=26).max()+df_orig.high.rolling(window=26).min())/2
    xdate = [x.date() for x in df_orig.index]
    
    plt.figure(figsize=(15,5))    
    plt.grid()
    plt.plot(xdate, df_orig.close,color="b",lw=1,linestyle="dotted",label="Close")
    plt.plot(xdate, df_orig.tenkan,label="Conversion line")
    plt.plot(xdate, df_orig.base,label="Base line")
    senkou1 = ((df_orig.tenkan+df_orig.base)/2).iloc[:-26]
    senkou2 = ((df_orig.high.rolling(window=52).max()+df_orig.high.rolling(window=52).min())/2).iloc[:-26]
    plt.fill_between(xdate[26:], senkou1, senkou2, color="blue",alpha=0.2,label="Cloud")
    plt.legend(loc=4)
    
    plt.xlim(xdate[0],xdate[-1])
    

    まとめ

    けっこういろいろな指標が考えられているのですね。 株価だけじゃなくて、他の時系列データへのヒントにもならないかと思い、調べましたが、とてもいい勉強になりました。

    ちょっと一つ移動平均線で試したいことがあるので、明日はそれをしようかなと思います。テーマは『ベイズ最適化』 それでは明日もお楽しみに!

    【Day-13】『Prophet入門』簡単に高精度を実現するFacebook謹製の時系列予測ライブラリ

    f:id:imslotter:20171213111429p:plain

    データ分析ガチ勉強アドベントカレンダー 13日目。

    仮想通貨がはやり始めて、チャートを見るようになった人も多いのではないでしょうか?

    チャートから予測をしたい

    という思いを持ちつつも、結構ハードルの高いのが時系列予測。 それをできるだけ簡単にできるツールがProphet

    自分の持っているドメイン知識を導入しながら、簡単に時系列データ予測を行うことができます。

    prophetとは

    資料

    • Completely automatic forecasting techniques can be brittle and they are often too inflexible to incorporate useful assumptions or heuristics.
    • Analysts who can produce high quality forecasts are quite rare because forecasting is a specialized data science skill requiring substantial experience.
    • 完全な予測は柔軟さに欠ける
    • 時系列予測は特別なスキルがいる

    から、良い予測が全然スケールしない。 なので、簡潔な記述必要な柔軟性を残し作ったのがProphet

    設計思想

    公式の設計思想がこちら

    • 簡潔な記述 : 正しい形で入力をすると、デフォルトで良い予測ができるようになっている
    • 柔軟性 : ユーザーの持つドメイン知識をいれることで、精度向上を測れる

    必要のないものは裏でやっちゃって、必要な部分だけ自由度を残すという感じですね。

    Prophetでできないこと

    これだけ言えば超絶便利だと思われるが、もちろんできないこともある。

    (認識間違っていたら教えてください。)

    • いろいろなアルゴリズムを試して比べるようなものではない*2
    • 予測のためのツール。時系列クラスタリングなどはできない。
    • 他変数を用いた予測(issueにあげられているので、そのうち対応されるかも。されてほしい)

    根底のアルゴリズム

    additive regression modelを利用している。具体的には下記のアルゴリズムを組み合わせている。

    • トレンドをlinearかlogistic曲線で近似する。また、Prophetではトレンドの変化点も自動で検知するようになっている
    • Fourier級数によって、年周期を抽出している
    • ダミー変数をもちいて週周期を抽出している
    • 重要な休日やイベントのリストをわれわれユーザー側で与えられる(ドメイン知識)

    インストール

    切り分けのため、別々にインストールするといいかも。*3

    # numpyとmatplotlibを最新版にしておくといい
    pip install -U numpy matplotlib
    # mcmc等サンプリング用ライブラリ
    pip install pystan
    # fbprophet
    pip install fbprophetc

    Prophet tutorial

    型の変換

    • pythonAPIでは、sklearnのモデルを踏襲しているため、fitpredictで学習と予測を行うことができる。
    • inputは必ず下記の型で行う
    変数 用途
    ds datetime 時系列情報(タイムスタンプ)
    y 数値 そのときの値
    %matplotlib inline
    import matplotlib.pyplot as plt
    import pandas as pd
    import numpy as np
    from datetime import datetime
    from fbprophet import Prophet
    df = pd.read_csv("AirportPassengers.csv",delimiter=";").dropna()
    def str2time(string):
        #datetimeに変換
        return datetime.strptime(string, '%Y-%m')
    # 読み込んだデータをdsとyに変える
    df["ds"] = df.Month.apply(str2time)
    df = df.drop(["Month"],axis=1)
    df["y"] = np.log(df.Passengers) #logをとる
    df.head()
    

    出力。dsyの形で持っている

       Passengers         ds         y
    0       112.0 1949-01-01  4.718499
    1       118.0 1949-02-01  4.770685
    2       132.0 1949-03-01  4.882802
    3       129.0 1949-04-01  4.859812
    4       121.0 1949-05-01  4.795791

    将来の予測

    下記の手順で、めちゃくちゃ簡単に予測までできる。面倒なパラメータ調整なんていらない。めちゃくちゃ簡単。

    • m = Prophet()で初期化(ドメイン知識はここで入れる!)
    • m.fit(df)で学習
    • make_future_dataframeで、dsのデータフレームを追加
      • freq : 予測する間隔(freq="M"ならMonthly。詳しくは公式 Sub-daily data )
      • periods : 何個dataframeを作るか
    • m.predct(future)で予測 : 出力はdataframe型
    • m.plot()でかなりオシャレなプロット

    今回は、下記要領で予測する

    • 1949年1月から1960年12年までのデータを入力にする
    • 1961年1月から5年間を予測する
    m = Prophet()
    # prophetに入れるには、dsとyのペアが必要
    m.fit(df)
    future = m.make_future_dataframe(periods=60,freq="M")
    forecast = m.predict(future)
    forecast[["ds","trend","yhat"]].tail()
    m.plot(forecast);
    plt.legend()
    

    出力がコレ

    かなりきれいな予測が出来ていそう。また、予測の誤差見積もりまでplotされている。色使いがFacebookっぽくてこれまたオシャレ感がある。

    Components

    アルゴリズムで紹介したように、トレンドや周期性についての分析結果は、m.plot_components()にて見ることができる。

    m.plot_components(forecast)
    

    特別なイベントなど

    ドメイン知識を入力し、精度向上を図ることができる。↓が参考になる。

    https://facebook.github.io/prophet/docs/seasonality_and_holiday_effects.htmlfacebook.github.io

    help

    その他使い方はhelp(Prophet)と打つ。

    まとめ

    一次元のデータを予測するなら、本当に簡単に分析が行えるので、オススメ。 株価でも突っ込んで分析をしてみようかな。

    株価といえば、明日はテクニカル分析で使われる指標について調べて実装してみます。オレオレ分析ツールを作りたいと思っているので。ではでは、また明日!

    *1:ドメイン知識の入れ方や、時系列予測に対する評価方法も書かれていて、良かった。

    *2:設計思想上その必要はない気がする

    *3:単にpip install fbprophetのみでもOKではありそう。

    【Day-12】時系列分析の良リソースまとめ&基礎チュートリアル

    データ分析ガチ勉強アドベントカレンダー 12日目。

    今までは、時間に依存しないデータについて取り扱ってきました。 しかし、世の中のデータは時間に依存したデータも多いのが事実です。

    時間に依存しないデータは、その分各データを独立に扱うことができますが、時系列データはそういうわけにはいきません。なので、なかなか難しいのです。

    今日は時系列のさわりをまとめて、また、時系列予測のチュートリアルをしていきます。

    参考にできるサイト

    時系列を学ぶにあたって、参考にした&参考にしたいサイトを列挙しておきます。 どれも充実した内容なので、一度は目を通しておくといいかもしれません。

    メタ的な記事

    pythonの実装も含めて

    ディープラーニング

    読んでる本、読みたい本を列挙しておく

    内容
    時系列解析入門 クラシカルな時系列データ分析が載っている本。数式をベースにきちんとした説明なのでよい。状態空間モデルのあたりは結構難しいが、たどり着きたいところ
    経済・ファイナンスデータの計量時系列分析 (統計ライブラリー) 未読。いいらしい。
    岩波データサイエンス Vol.6 時系列に関する最新トピックを集めた本。かなり面白いが、それなりの素養が必要
    詳解 ディープラーニング ~TensorFlow・Kerasによる時系列データ処理~ KerasとTensorflowを使ってRNNを実装する本。LSTMの説明などもされていて、ディープラーニングのモデルを実装しつつ勉強ができる。

    読みたい論文(積読)

    (見出しを押せば論文にとびます)

    Toeplitz Inverse Covariance-Based Clustering of Multivariate Time Series Data

    2017年KDDのBestPaperを受賞している。時系列を分析して状態をクラスタリングするという旨の論文。

    www.youtube.com

    WAVENET: A GENERATIVE MODEL FOR RAW AUDIO

    CNNで時系列データ分析を実現した論文。Google音声合成などにも使われていて、高い性能を実現している

    日本語の輪読スライドも、ぜひとも参考にしたいところ。

    www.slideshare.net

    基礎チュートリアル

    上記で紹介した記事A comprehensive beginner’s guide to create a Time Series Forecast (with Codes in Python) を参考にしながら、実際に時系列の予測をしてみる。

    ところどころライブラリのバージョンが古かったりしたので、修正&大事そうなところは加筆した。内容は↓

    • データの用意
    • 定常性チェック
    • 定常状態のデータに変換
    • 時系列の予測

    いつものようにgithubにあげている(day12)

    github.com

    データの用意

    pandasを用いてデータを用意する 使うデータは International Airport Passengers

    定常性のチェック

    • 平均や分散がずっと一定ならstarionary(定常)といえる
    • 多くの時系列モデルは定常性を仮定しているので、定常性のチェックは重要。

    使用しているデータを見ると、増加傾向にあることがわかる(上図)。これだと平均や分散が一定とならないように見える。 ちゃんと調べるために、下記の処理を施してみる

    1. moving averageやmoving distanceをとってみる(sliding window)今回は、季節変動的なのも吸収したいから12か月分
    2. Dickkey Fuller Testを用いる。定常性の統計検定である。詳細

    移動平均や移動標準偏差には、pandasの関数rollingを用いる。(pandas rolling) windowは窓幅で、window=12はこのデータでは一年分の平均となる。これで、簡単に移動平均をとることができる。

    rolmean = timeseries.rolling(window=12).mean()
    

    移動平均/標準偏差をとった結果は↓

    DF検定には、statsmodelsという時系列用統計ライブラリを使う。statsmodelsは時系列を取り扱うときに何かと便利なので、扱えるようにしておきたい(statsmodels)

    from statsmodels.tsa.stattools import adfuller

    で、検定をしてくれる。検定結果は↓

    Results of Dickey-Fuller Test:
    Test Statistic                   0.815369
    p-value                          0.991880
    #Lags Used                      13.000000
    Number of Observations Used    130.000000
    Critical Value (1%)             -3.481682
    Critical Value (5%)             -2.884042
    Critical Value (10%)            -2.578770

    この表の見方だが、Test StatisticとCritical Valueを見比べる Test Statisticの値が、Critical Valueを下回らないと定常とはいえない。

    たとえば、Test Statisticが-3.0だったとすると、95%以上の信頼を持って定常といえると結論付けることができる。ともかく、このままだとこいつはまったく持って定常ではないデータなのである。

    時系列データを定常状態にしたい

    変動要因にはいろいろある

    • トレンド
    • 季節変動
    • ノイズ
    • 小変動(周期性あり)

    定常性に響くトレンドをうまく取り除き、定常状態にしてやりたい。 下記の要領で行っていく

    1. positiveトレンドのときは、大きい値にペナルティを与えるような変換をかける。さまざまな変換がある。

      • logをとる
      • logit変換{$z_n = \log{\frac{y_n}{1-y_n}}$}
      • 差分
      • ルート、三重根をとる
    2. トレンド成分を取り出す

      • Aggregation : 週/月ごとにまとめた平均を取る
      • Rolling : windowをスライドさせながら平均を取る(移動平均)
      • Polynomial Fitting : 回帰でフィッティングさせる
    3. Originalからトレンド成分を引く
    4. テストをしてみる。

    の要領でやっていく。今回はlogをとって、移動平均にかけて、Originalから引いてみた。 移動平均には、普通の平均をとるのと、最近の時間ほど重みを強くかけるexpweighted moving average**の二種類をとった。後者は時間幅を厳格に設定しなくてよいため、クリティカルな窓幅が想像できないときなどに、柔軟に扱える。

    moving averageによってトレンドを求めた結果がこちら

    このトレンド成分(weighted moving average)をOriginalから引いたのがこちら

    この結果を検定にかけると、下記の結果になる

    Results of Dickey-Fuller Test:
    Test Statistic                  -3.601262
    p-value                          0.005737
    #Lags Used                      13.000000
    Number of Observations Used    130.000000
    Critical Value (1%)             -3.481682
    Critical Value (5%)             -2.884042
    Critical Value (10%)            -2.578770

    99%以上定常状態といえる。このようにして定常性を消すと良い。

    時系列の予測

    最後に予測。 予測によく使われるARIMA(Auto-Regressive Integrated Moving Average)の説明をしながら、実際に使ってみる。 ARIMAは、ARモデルと、MAモデルを統合したモデルとなっている。

    ARIMAには、下記のパラメータを用いる

    • p (number of AR terms) : 自己回帰の窓幅を何個前まで取るか
    • q (number of moving average) : 移動平均の窓幅を何個前まで取るか
    • d (number of differences) : 一次の差分か二次の差分か

    p=0だとMA、q=0だとARモデルとできる。

    pとqを決めるための大事な量がACF(Autoorrelation FunctionPACF(Partial ACF)

    • p : ACFのupper confidential intervalがクロスするラインを選ぶ
    • q : PACFupper confidential intervalがクロスするラインを選ぶ

    と良いらしい。これもstatsmodelsで利用可能。

    from statsmodels.tsa.stattools import acf, pacf
    

    *1 よさげなpとqがわかったら、ARIMAにかけてみる。 これもstatsmodelsで利用ができる。

    • モデル初期化はARIMA(p,d,q)の順番
      • p=0でMAモデル
      • q=0でARモデル
    • sklearnのように、fitでモデル作成
    ts #時系列データ
    from statsmodels.tsa.arima_model import ARIMA
    model = ARIMA(ts, order=(2, 1, 0))  
    model.fit()
    

    いろいろとわちゃわちゃしながら*2の結果がこちら

    定量評価には、RMSE(Root Mean Square Error)を使ってみた。自分で簡単に実装できるが、復習もかねてsklearn.metricsを使ってみる。

    from sklearn.metrics import mean_squared_error
    print("ARIMA:", np.sqrt(mean_squared_error(df.Passengers, df.ARIMA)))
    print("MA", np.sqrt(mean_squared_error(df.Passengers, df.MA)))
    print("AR", np.sqrt(mean_squared_error(df.Passengers, df.AR)))
    

    今回はARが最も良かった模様?

    ARIMA: 50.2956404427
    MA 43.9652476837
    AR 20.4183022135

    まとめ

    時系列分析事始ということで、調べたいリソースを集約するのと、はじめのチュートリアルまでを行った。 statsmodelsをとりあえず知っておけば、ぱぱっとした分析には良いなという感じなので、習得しておきたいところ。

    しかし最近、facebookが時系列ツールを出したので、そちらのキャッチアップもしておきたい。 なので、明日はfacebook謹製の時系列予測ツールprophetについて触れたい。

    明日もお楽しみに!

    *1:でも結局pとqは試行錯誤で選ぶ気もしている

    *2:詳しくはgithbをご覧ください

    【Day-11】機械学習のチートシートを眺めたり、比べてみたり

    f:id:imslotter:20171210225755p:plain

    データ分析ガチ勉強アドベントカレンダー 11日目。

    モデルを選び、試行錯誤しながら作っていく、そんな過程まで勉強してきました。 実装寄りの内容になったので、ここで一度機械学習界隈を俯瞰してみようと思いまして、調べると出てくる有名なチートシート

    について、ここでまとめておきます。

    用途ベースのMicrosoft Azure Machine Learning Cheat Sheet

    Microsoftが提供しているチートシート用途別にできるだけ分岐を行って、その後でデータの性質別に分けている。

    用途 アルゴリズム
    データの隠れた構造を見つけたい クラスタリング
    カテゴリを予測したい 分類
    値を予測したい 回帰
    普通じゃないデータを予測したい 異常検知

    データベースのMachine Learning Cheat Sheet (for scikit-learn)

    こちらはpython機械学習を行うときには必須のツールsklearnチートシート

    どういうデータかというところに着目してチートシートが描かれている。

    その中で、データの種類によって4種類に分けている。ニュアンスは下表な感じ

    データの種類 アルゴリズム
    カテゴリのラベルを持っている 分類
    カテゴリラベルを持っていない クラスタリング
    値のラベルを持っている 回帰
    どうしようもないデータを見やすくしたい 次元削減

    両者の違い

    上の説明を見ればわかるが、Microsoftチートシートは用途ベースで話を進めている。なので、どういうことをしたいかって言うのがわかりやすい。また、ほとんど定性的な説明で書かれているので、イメージがしやすい。

    一方でsklearnはデータの中身を見ながら分岐している。チェーンソーで紙を切るような無駄を起さないために、また良い妻で立っても計算が終わらないなどの悲劇を起さないために、なるべく簡単な学習器から初めてうまくいかなかったら高度な学習器をっていう感じになっている。

    ぱっと見タスクは何かを把握するにはMicrosoftチートシートだが、データ分析従事者にとってはsklearnのチートシートを見るのがいいのかなっていう気がした。各アルゴリズムの中身を見ていると

    • (回帰)一部の特徴量が効いているっぽい : LASSOなどによるスパース正則化を導入
    • (分類)テキストデータ : ナイーブベイズ
    • サンプルデータが大きい非線形分類 : カーネル近似

    など、割と細かく分かれている。これをベースにアルゴリズムを勉強して、「なぜそういう分岐になっているのか」を理解していくのもよさそう。

    見て判断できるようにしておく

    上記チートシートアルゴリズムを選ぶ際の参考になる。コレに加えて、自分で実験しながら視覚的&数値的にその学習器の性質を判断できると良いのかなと思う。Day-9で少し触れたが、評価したり可視化したりできていると、たとえば分類ならどんな分類面の引き方をするのかな~がわかる。

    www.procrasist.com

    いつものGithubに、データを投げるといろんな学習器を調べるスクリプトを作ってみた。

    github.com

    まだ改良中だが、これからもう少し整備したい。

    分類

    今回は簡単な分類問題と、三日月形の分類で、その識別境界を色付けしてみた。濃いほど自信があるように色をつけてある。

    blobs

    線形分離がほぼできる分布 LogisticRegressionやLinearSVMなどのシンプルな分類器が一番わかりやすそう。 一方で、RBFカーネルはちょっとぐにゃぐにゃとした分離選を作りすぎかなというようにも感じられる。 加えて、アンサンブル学習系(Decision Tree/ Adaboost/ RandomForest)なども斜めに選を引かずかくかくとした線になっていて少し苦しそう*1

    三日月形

    三日月形になると、線形分類器たちは太刀打ち不能になってしまう。一方で、GaussianProcessなどの非線形分類器は、かなりきれいな分離面を引いているのがわかる。ただGaussianProcess等は計算コストが大きいので、そのあたりも考慮しながらの運用になりそう。

    クラスタリング

    クラスタリングも、パラメータの有無だったり、クラスタ形成の仕方だったりで、いろいろと結果に違いが生じる

    blobs

    簡単なクラスタリング問題では、ほとんどのアルゴリズムが正しく機能している。 一方で、DBSCANだけはクラスタが多くなってしまっているようだ。

    三日月形

    じゃあDBSCANが悪いのかといわれると、実はそんなことはない。非線形な分布クラスタリングでは威力を発揮することも多い。下図だと、DBSCANだけが狙ったクラスタリングを達成できている。 閾値の調整が結構めんどくさいが、パラメータ次第では協力に機能する。

    まとめ

    きょうはチートシートを眺めたり、アルゴリズムを俯瞰したりしてみた。 結局時と場合によることが多いので、そのときにより良い判断ができるような基準を自分の中でもっておくことが必要かなと思う。

    などを総合的に判断できるようになると、データサイエンティストに近づけるのかなと思う。

    結構長い間sklearnメインで話を進めてきたけど、データ分析の基礎ができたところで、明日からは少し別の話題に触れようと思う。明日は時系列データについて書いていきたい。それではまた明日。

    *1:座標変換すれば何とかなるのかな?

    【Day-10】Cross Validationとパラメータサーチでモデルの調整

    f:id:imslotter:20171210183652p:plain

    データ分析ガチ勉強アドベントカレンダー 10日目。

    • データを集め、前処理を行い、学習をする。
    • どういう学習器が良いのかの評価基準

    の勉強までできた。でも、データがあって、評価基準がわかっていても、どうやって評価すればいいかについてはまだあまり触れていない。

    もうちょっと具体的に言うと、データと学習器があって、どういう風に教師データとテストデータを分ければいいのか。この方法について、実はまだ何も言っていない。

    そこで登場するのが、sklearn.model_selection

    いよいよ最後の仕上げ。良い学習器、良いパラメータの探り方について学ぶ

    Validationとは

    モデルを作ってデータを分析する際、データをtrain data, validation data, test dataに分ける。 最終的なテストをする前に、訓練データの中を分割してテストを回すことで、パラメータ調整を行うために用いられる。

    sklearnでは、originalデータをtrainとtestに分けるところをsklearn.model_selection.test_train_splitで行い、Validationデータに分けるところにもsklearn.model_selection.cross_val_scoreなどが用意されている。

    また、それにしたがって学習器のパタメータを変えながら評価していく仕組みも用意されていて、sklearn.model_selection.GridSearchCVで実現される。(後述)

    skleanのCross Validation

    cross_val_score

    sklearnで最も簡単にCross Validationをするには、cross_val_scoreという関数を用いる。

    `cross_val_score(clf, data, target, cv=5, scoring="accuracy")`
    
    変数 役割
    clf 学習器
    data データ
    target ラベル
    cv クロスバリデーションの回数
    scoring 評価する指標(参照:metrics)

    scoringには、下記の値を用いることができる。

    ['accuracy', 'adjusted_rand_score', 'average_precision', 
    'f1', 'f1_macro', 'f1_micro', 'f1_samples', 'f1_weighted', 'neg_log_loss', 
    'neg_mean_absolute_error', 'neg_mean_squared_error', 'neg_median_absolute_error', 
    'precision', 'precision_macro', 'precision_micro', 'precision_samples', 'precision_weighted', 
    'r2', 'recall', 'recall_macro', 'recall_micro', 'recall_samples', 'recall_weighted', 
    'roc_auc']
    

    各評価値が何なのかを詳しく知りたい方は、Day9を参照

    出力は、cvの回数分の、評価値の配列となっている。例として、digitsデータをRandomForestで取り扱ってみる。

    from sklearn.datasets import load_digits
    from sklearn.ensemble import RandomForestClassifier
    from sklearn.model_selection import cross_val_score
    digits = load_digits()
    clf = RandomForestClassifier()
    scoring = "f1_macro"
    scores = cross_val_score(clf, digits.data, digits.target, cv=5, scoring=scoring)
    print("{}:{:.3f}+/-{:.3f}".format(scoring, scores.mean(), scores.std()))
    

    結果

    f1_macro:0.901+/-0.031

    pipline

    cross validationを行う際に、sklearn.pipelineを使うと、処理が簡潔に書ける。(参考)

    例 : 正規化した後の効果をcross_validationしたい

    # 通常
    from sklearn import preprocessing
    X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, train_size=0.8)
    scaler = preprocessing.StandardScaler().fit(X_train)
    X_train_transformed = scaler.transform(X_train)
    clf = SVC(C=1).fit(X_train_transformed, y_train)
    X_test_transformed = scaler.transform(X_test)
    print(clf.score(X_test_transformed, y_test))
    # pipeline
    from sklearn.pipeline import make_pipeline
    # 処理を流す順に関数を書いていく
    scoring = "f1_macro"
    clf = make_pipeline(preprocessing.StandardScaler(), SVC(kernel="linear",C=1))
    scores = cross_val_score(clf, digits.data, digits.target, cv=5, scoring=scoring)
    print("{}:{:.3f}+/-{:.3f}".format(scoring, scores.mean(), scores.std()))
    

    pipelineを使うと、複数処理を一括でCVにかけられるので、とてもよい。

    cross_validate 関数

    sklearn 0.19.1から、cross_validateが用意された。

    さっきまでのがcross_val_scores関数。ややこしいが少し違う

    • 評価に複数の指標を考慮できる
    • テストスコアに加えて、学習の時のスコア、学習時間、テストの時間などを算出してくれる。

    つまり、より強力なcross validationを行える。

    from sklearn.model_selection import cross_validate
    from sklearn.metrics import recall_score
    from sklearn.svm import SVC
    scoring = ["f1_macro", "recall_macro"]
    clf = SVC(C=1)
    scores = cross_validate(clf, digits.data, digits.target, scoring=scoring, cv=5)
    for key,value in scores.items():
        print("{}:{:.2f}+/-{:.2f}".format(key, value.mean(), value.std()))
    

    結果はこちら

    fit_time:0.51+/-0.01
    score_time:0.12+/-0.00
    test_f1_macro:0.50+/-0.04
    train_f1_macro:1.00+/-0.00
    test_recall_macro:0.45+/-0.04
    train_recall_macro:1.00+/-0.00

    Validation iterators

    CrossValidationを行うときの、データセットの分け方にも、いろんなバリエーションがある。 sklearn.model_selection内には、いろんなCVの分け方の方法が載っている。

    • i.i.d data
      • K-fold
      • Repeated K-Fold
      • Leave One Out (LOO)
      • Leave P out (LPO)
      • Shuffle & Split
    • iterators with stratification based on class labels(サンプリングとか)
      • Stratified k-fold
      • Stratified Shuffle Split
    • Grouped data i.i.dデータと基本おなじ

    代表的なのが、K-Fold法とLeave One Out法。両者を図にするとこんな感じ

    今回は、K-Fold法でのCrossValidationを実装してみる

    import numpy as np
    from sklearn.model_selection import KFold
    from sklearn.datasets import load_digits
    from sklearn.metrics import accuracy_score
    # K-Fold交差検定
    iterater = KFold(n_splits=5)
    results = []
    for train_indexes, test_indexes in iterater.split(digits.data):
    #     print(train_indexes, test_indexes)
        X = digits.data[train_indexes]
        y = digits.target[train_indexes]
        clf = RandomForestClassifier()
        clf.fit(X,y)
        pred_y = clf.predict(digits.data[test_indexes])
        ac = accuracy_score(digits.target[test_indexes], pred_y)
        results.append(ac)
    results = np.array(results)
    print("KFold accuracy: {:.2f}+/-{:.2f}".format(results.mean(),results.std()))
    

    出力

    KFold accuracy: 0.90+/-0.04

    パラメータを調整する

    機械学習は複雑なので、モデルの持つパラメータをどのように調節すべきか難しい。

    sklearnでは、下記の方法でパラメータ探索が行えるようになっている(参考 : Tuning the hyper-parameters of an estimator)

    • グリッドサーチ : パラメータを総当り式で調べる
    • ランダム選択 : ランダムにパラメータをサンプリングして試す
    • モデルに合わせた賢いCV*1

    他に、変数最適化の問題として、ある指標を評価関数にして解くというのもよさそう。Gaussian Process Optimizationなどがよさげ (参考)

    今回は、グリッドサーチでの方法を載せる

    パラメータ

    SVCのパラメータである、C, kernel, その他変数について調整する。調整する パラメータは下記のように、キーとリストで保持する。

    param_grid = [
        {'C': [1, 10, 100, 1000], 'kernel': ['linear']},
        {'C': [1, 10, 100, 1000], 'kernel': ['rbf'], 'gamma': [0.001, 0.0001]},
        {'C': [1, 10, 100, 1000], 'kernel': ['poly'], 'degree': [2, 3, 4], 'gamma': [0.001, 0.0001]},
        {'C': [1, 10, 100, 1000], 'kernel': ['sigmoid'], 'gamma': [0.001, 0.0001]}
        ]
    

    関数にはsklearn.model_selection.GridSearchCVを用いる。引数は

    • 学習器のインスタンス(今回はSVC)
    • パラメータグリッド(上記param_grid)
    • cv : cross validationの回数
    • scoring : 評価する指標

    実装

    実装は下記の通り

    from sklearn.model_selection import GridSearchCV
    import pandas as pd
    clf.get_params()
    svc = SVC()
    scoring = "f1_macro"
    param_grid = [
        {'C': [1, 10, 100, 1000], 'kernel': ['linear']},
        {'C': [1, 10, 100, 1000], 'kernel': ['rbf'], 'gamma': [0.001, 0.0001]},
        {'C': [1, 10, 100, 1000], 'kernel': ['poly'], 'degree': [2, 3, 4], 'gamma': [0.001, 0.0001]},
        {'C': [1, 10, 100, 1000], 'kernel': ['sigmoid'], 'gamma': [0.001, 0.0001]}
        ]
    clf = GridSearchCV(svc, param_grid,cv=4)
    clf.fit(digits.data, digits.target)
    df = pd.DataFrame(clf.cv_results_)
    df_scored = df.sort_values(by=["rank_test_score"])[["params","mean_test_score","std_test_score","mean_fit_time"]]
    

    結果を見る

    結果は、pandasのデータフレームに保管可能になっていて、簡単に好きな値でソートができる。今回はテストスコアの良い順にソートを行ってみた。(pandasの扱いが苦手な人はDay4をご覧ください)

    こんな風に、どういうパラメータがどういう結果だったかが保存されているので、より良いパラメータを選ぶことができる。

    まとめ

    今日は学習器の調整に重要なCross Validationとパラメータサーチの方法について触れた。 実際の業務の場では、パラメータを調整する機会もとても多いので、知っておくと良い。また、少し述べたが、もう少し効率的なパラメータ調整方法についても勉強しておきたい。

    機械学習の一連の流れがコレでざっと洗えたと思う。明日はscikit-learnのチートシートについてでも触れようかな。明日もお楽しみに!

    *1:ドキュメントによればあるらしい...?

    【Day-9】機械学習で使う指標まとめ(実装編)

    f:id:imslotter:20171209195625p:plain データ分析ガチ勉強アドベントカレンダー 9日目。

    データを学習器に入れるところまではできた。後は学習させるだけ! だが、学習器といってもたくさんある。どういう学習器を選べばよいのだろうか。

    そのためにはモデルをうまく評価するしくみを作らなければならない。 今回は、sklearn.metricsを見ながら、学習の評価をどのように行うかについて学ぶ

    機械学習に使う指標総まとめ(教師あり学習編)

    学習器の評価には、さまざまな指標が用意されている。 以前、教師あり学習については全力でまとめている。なので、そちらを見ていただければと思う。

    www.procrasist.com

    しかし、こちらの記事では、数式も多く、実感もわかないだろう。なので今回はこの中の指標の実装も行ってみる。

    sklearn.metricsによる定量指標

    sklearn.metricsには、たくさんの指標が用意されている

    分類の定量指標

    先ほどの分布を分類した結果を評価したい。参考にしたのはこちらの公式ページ。 用意するものは、下記のものである。カッコ内は私が実装で用いた変数

    • テストデータの正解ラベル(y_test )
    • 学習器が判定した結果(pred_y )
    • そのときの判定スコア(pred_y_score)※

    ※判定スコアは、必要な指標と必要でない指標がある。使わない指標は、ただ分類結果だけをみて判定するのでわかりやすい。一方で判定スコアを使う場合は、指標がすこしわかりにくくなる一方で、いわゆる分類器の自信を考慮しながらの算出となるので、より詳細に分類器を評価している。

    ここで、sklearn.metricsに入っている主要な分類指標を見てみる

    クラス 分類指標 判定スコアの
    使用
    備考
    confusion matrix -- なし 実際と予測のずれ具合を見る
    accuracy Accuracy なし どのくらいあたったか
    classification_report f値, recall, precision なし よく使う指標をまとめてくれてて良い
    hamming_loss hamming loss なし f:id:imslotter:20171209182300p:plain
    jaccard_similarity jaccard類似度 なし f:id:imslotter:20171209182351p:plain
    log_loss logarithm loss あり 使い方あんまり自信がない公式
    precision_recall_curve precision-recall曲線 あり トレードオフ具合をみる
    roc_auc ROC曲線 あり 閾値判定とか

    よく使いそうだなと思ったのが、accuracy. classification_report. roc_auc。主要な指標がきれいにまとまっている

    sklearn.datasets.load_breast_cancer()の学習を簡易的に行い、上記の3つの指標の入力の仕方と出力の仕方を学ぶ。

    コードをみる
    import matplotlib.pyplot as plt
    from sklearn.datasets import load_breast_cancer
    from sklearn.linear_model import LogisticRegression
    from sklearn.model_selection import train_test_split
    from sklearn import metrics
    breast_cancer = load_breast_cancer()
    clf = LogisticRegression()
    X_train, X_test, y_train, y_test = train_test_split(breast_cancer.data, breast_cancer.target,train_size=0.8)
    clf.fit(X_train,y_train)
    # 結果を返す
    pred_y = clf.predict(X_test)
    # スコアを出すための関数。
    # ref : [decision_functionとの違いは?](https://stackoverflow.com/questions/36543137/whats-the-difference-between-predict-proba-and-decision-function-in-sklearn-py)
    pred_y_score = clf.predict_proba(X_test)[:,1]
    print("accuracy\n", metrics.accuracy_score(y_test, pred_y))
    print("classification report\n",metrics.classification_report(y_test,pred_y))
    # Compute ROC curve and ROC area for each class
    fpr, tpr, _ = metrics.roc_curve(y_test, pred_y_score,pos_label=1)
    roc_auc = metrics.auc(fpr, tpr)
    print("auc area", roc_auc)
    # ROC
    plt.figure(figsize=(5,5))
    plt.plot(fpr, tpr, color='darkorange', lw=2, label='ROC curve')
    plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('ROC curve: \nAUC={:.3f}'.format(roc_auc))
    plt.show()
    

    結果はこんな感じ

    f:id:imslotter:20171209184215p:plain

    回帰の定量指標

    正解が数値で与えられる回帰にも、いくつか指標がある。これもsklearn.metrics内で実装されている。

    sklearn.datasets,load_bostonデータセットを使って、試してみる

    コードをみる
    from sklearn.datasets import load_boston
    from sklearn.linear_model import LinearRegression
    from sklearn import model_selection, metrics
    import numpy as np
    boston = load_boston()
    reg = LinearRegression()
    X_train, X_test, y_train, y_test = model_selection.train_test_split(boston.data, boston.target,train_size=0.8)
    reg.fit(X_train,y_train)
    pred_y = reg.predict(X_test)
    print("Variance score", metrics.explained_variance_score(y_test, pred_y))
    print("MAE : Mean absolute error", metrics.mean_absolute_error(y_test, pred_y))
    print("MSE : Mean squared error", metrics.mean_squared_error(y_test, pred_y))
    print("RMSE: Root mean squared error", np.sqrt(metrics.mean_squared_error(y_test, pred_y)))
    print("R2 score", metrics.r2_score(y_test, pred_y))
    

    出力は

    Variance score 0.637022922972
    MAE : Mean absolute error 3.16480595297
    MSE : Mean squared error 19.7557043856
    RMSE: Root mean squared error 4.44473895584
    R2 score 0.633791811542

    クラスタリング定量指標

    Evaluating the performance of a clustering algorithm is not as trivial as counting the number of errors or the precision and recall of a supervised classification algorithm. (clustering performance evaluation)

    とあるように、クラスタリングこそ絶対的な定量指標が存在しない。分け方なので、用途や意味によって変わってくるから。 けれど、一応指針みたいなのはあって、クラスタ内の散らばりの少なさクラスタ間の距離の長さをうまく定式化することで、定量化する。

    上記サイトにはクラスタリングの指標も多数実装されている。*1

    デモ実装(分類学習器ごとの差)

    最後にデモを作る。

    この図は、三日月形に分布するデータを分類しようとした結果を分類器別に分けてみたものである。 *2

    線形、非線形、アンサンブル学習など、学習器の性質によってかなり異なることが見て取れて面白い。

    こうやって分離状況を可視化していくと、大体どういう雰囲気で学習器が分布しているかはわかる。 でも、それでは定量的な指標とはいえない。ので、先ほどまでの指標を見て、定量化してみた

    詳細は例のごとくgithubにあげている。

    github.com

    結果1(表)

    大きく結果の順位が変わるようなことはないが、精度だけじゃなく、precisionとrecallなどを見ると開きがあるもの、ないものなど違いがあったりして、分類器の個性を詳細に知ることができる。

    Algorithm accuracy precision recall f1-score ham. loss jaccard sim AUC
    LogisticRegression 0.83 0.82 0.83 0.82 0.17 0.83 0.92
    Nearest Neighbors 0.95 0.91 1 0.95 0.05 0.95 0.95
    Linear SVM 0.8 0.8 0.77 0.79 0.2 0.8 0.92
    RBF SVM 0.94 0.89 1 0.94 0.06 0.94 0.98
    Gaussian Process 0.95 0.91 1 0.95 0.05 0.95 0.98
    Decision Tree 0.92 0.92 0.92 0.92 0.08 0.92 0.92
    Random Forest 0.89 0.88 0.9 0.89 0.11 0.89 0.93
    Neural Net 0.83 0.82 0.83 0.82 0.17 0.83 0.92
    AdaBoost 0.86 0.84 0.88 0.86 0.14 0.86 0.92
    Naive Bayes 0.82 0.81 0.81 0.81 0.18 0.82 0.92
    QDA 0.83 0.82 0.83 0.82 0.17 0.83 0.92

    結果2(図)

    代表的な3種の分類器でROC曲線とPrecision-Recall曲線を引いた。 両者とも、トレードオフのグラフになっているので、覆われた領域の面積が大きいほど、いい分類器であるといえる。 今回のデータに対しては、非線形分類器に軍配が上がりそうだ。

    まとめ

    以前まとめていた指標、理論ではわかっていても実際に使って動きをみてみると、知識として定着する。 いろいろと指標はあるが、sklearn.metricsを用いれば簡単に実装が可能なので、勉強したい方は是非一度手を動かしてみてはいかがでしょうか。

    ちょっと長くなりましたが、今日はこの辺で。 明日はモデル選択(Cross Validation)とパラメータサーチ!

    *1:余裕があれば加筆します...w

    *2:Classifier Comparisonを一部改変してつくった

    PROCRASIST