雑記: AUC の話とその scikit-learn での計算手順の話

AUC というものを算出しろといわれることがあると思います。でも幸いなことに scikit-learn で算出できます。

例えば以下のような正解ラベルが付いたデータのそれぞれに対して、あるモデルが以下のようなスコアを出しているとします。正解ラベルを「スパムメールか否か」、モデルのスコアを「モデルが推定したスパムメールらしさ」などと考えてもいいです。

真の正解ラベルモデルの出力スコア
02
01
02
14
12
01
13
15
このモデルのスコアの AUC を得るには、上のデータを sklearn.metrics.roc_auc_score に渡すだけです。あるいは sklearn.metrics.roc_curve に渡してからその出力を sklearn.metrics.auc に渡しても同じです( AUC とはその名の通り「Area Under Curve=曲線の下の面積」なので、後者は一旦曲線を出していることになります )。

from sklearn import metrics

list_label = [0, 0, 0, 1, 1, 0, 1, 1]
list_score = [2, 1, 2, 4, 2, 1, 3, 5]

auc = metrics.roc_auc_score(list_label, list_score)
print(auc)

fpr, tpr, thresholds = metrics.roc_curve(list_label, list_score)
auc = metrics.auc(fpr, tpr)
print(auc)
0.9375
0.9375

AUC は 0.9375 となります。この AUC がどんな曲線の下の面積だったのかを知るには、metrics.roc_curve の出力をプロットしてみるとわかります。

%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 6, 4
rcParams['font.family'] = 'Ume Gothic O5'
rcParams['font.size'] = 16

print("fpr: ", fpr)
print("tpr: ", tpr)
print("thresholds: ", thresholds)

plt.plot(fpr, tpr, marker="o")
plt.xlabel("fpr")
plt.ylabel("tpr")
plt.show()
fpr:  [0.  0.  0.5 1. ]
tpr:  [0.25 0.75 1.   1.  ]
thresholds:  [5 3 2 1]

f:id:cookie-box:20190209172457p:plain

fpr、tpr なるものをプロットし、各点を直線で結ぶとこのような折れ線が描けます。そして、この折れ線の下の面積は確かに 0.9375 になっていることがわかります(※ 上図では補っていませんが、折れ線は本当は原点から始まっていると考えてください)。
しかし、fpr、tpr、thresholds というのは何なのでしょうか。それは thresholds というリストに入っている値のそれぞれを閾値としたときのモデルの判定を下のように書き出してみるとみえてきます(下の表では、8つのデータをスコア降順に並べ直しています)。
正解ラベルスコア▼判定(閾値=5)判定(閾値=3)判定(閾値=2)判定(閾値=1)
1:スパム5スパム(正)スパム(正)スパム(正)スパム(正)
1:スパム4非スパム(誤)スパム(正)スパム(正)スパム(正)
1:スパム3非スパム(誤)スパム(正)スパム(正)スパム(正)
0:非スパム2非スパム(正)非スパム(正)スパム(誤)スパム(誤)
0:非スパム2非スパム(正)非スパム(正)スパム(誤)スパム(誤)
1:スパム2非スパム(誤)非スパム(誤)スパム(正)スパム(正)
0:非スパム1非スパム(正)非スパム(正)非スパム(正)スパム(誤)
0:非スパム1非スパム(正)非スパム(正)非スパム(正)スパム(誤)
真のスパムメールのうち
スパムと判定されている割合
25%75%100%100%
真の非スパムメールのうち
スパムと判定されている割合
0%0%50%100%
上の表の下から2番目の行が tpr 、1番下の行が fpr に相当していることがわかります。しかし、この下の面積にどんな意味があるのでしょうか。
上の表の下から2番目の行は閾値をどんどん下げていくとどれだけスパムメールがカバーされていくか」と解釈できます。なので、0% から始まり、やがて 100% に至ります(本当は上の表の左端に 0% が補われます)。他方、一番下の行は閾値をどんどん下げていくとどれだけ非スパムメールまで巻き込んでカバーされてしまっていくか」と解釈できます。これも 0% から始まり 100% に至ります(やはり上の表の左端に 0% が補われます。今回のケースはたまたま補わなくても一番左が 0% ですが、いつもそうではありません)スパムメール判定モデルに求められることは、閾値をどんどん下げていくと、先にスパムメールだけがカバーされていく(非スパムメールを巻き込まずに)」、そんなスコアを出力することです。
  • 例えば完璧なモデルなら、閾値を下げていったときに先に完全にスパムメールをカバーし切って、その後非スパムメールがカバーされ始めます。つまり、fpr=0% のまま(非スパムメールを1件も巻き込まないまま)tpr=100% に到達するので(全てのスパムメールをカバーするので)、このときの fpr-tpr 曲線は (0%、0%)→(0%、100%)→(100%、100%)の形を描きます。この曲線の下の面積は 1 です。
  • 逆にスパムメールと非スパムメールを全然識別できないモデルなら、閾値を下げていったときにスパムメールも非スパムメールも同じ割合ずつカバーしていってしまいそうです。このときの fpr-tpr 曲線は (0%、0%)→(100%、100%)を直線に結びます。この曲線の下の面積は 0.5 です。
なので、AUC が 1 に近いスコアほどよいスコアであるといえそうです(AUC の値は非スパムメールよりスパムメールに大きいスコアが付いている期待値などと解釈もできます。参考: 昔自分が書いた非常に読みづらい記事(※)
AUC が計算できたのはいいんですが、時には scikit-learn の力を借りることができないこともあると思います。でも scikit-learn と値が合わないと気持ち悪いと思います。上の流れをみると自分でも適当に実装できそうな気がしますが雰囲気でやるのはよくないと思います。ちゃんとコードをみます。sklearn.metrics.roc_auc_score は以下です。

これをみると結局 sklearn.metrics.roc_curve と sklearn.metrics.auc がよばれていることがわかります。roc_curve は以下です。

この引数をみて1つ解決した疑問があります。上の例では roc_curve の戻り値の thresholds は [5 3 2 1] となっており 4 が含まれていませんでした。しかし 4 というスコアのデータもあるので、「閾値=4」での fpr、tpr も計算するべきではないかという気もします。それは roc_curve の drop_intermediate という引数がデフォルトで True だったために、ROCカーブの形状を変えない座標の fpr、tpr が省かれていたということがここにきてわかります。試しに、drop_intermediate=False としてもう一度実行してみます。

fpr, tpr, thresholds = metrics.roc_curve(list_label, list_score, drop_intermediate=False)
print("fpr: ", fpr)
print("tpr: ", tpr)
print("thresholds: ", thresholds)

plt.plot(fpr, tpr, marker="o")
plt.xlabel("fpr")
plt.ylabel("tpr")
plt.show()
fpr:  [0.  0.  0.  0.5 1. ]
tpr:  [0.25 0.5  0.75 1.   1.  ]
thresholds:  [5 4 3 2 1]

f:id:cookie-box:20190210162832p:plain

thresholds に 4 が含まれ、ROCカーブの座標が先ほどより1点増えました。しかし、この座標はあってもなくてもROCカーブの下の面積を出す上で影響がない点であったこともわかります。ROCカーブを軽くするためにデフォルトではこのような点は省かれているようです。
それで roc_curve ではまずサブ関数の _binary_clf_curve がよばれて、閾値をどんどん下げていったときに検出されるスパムメールと非スパムメールの個数を出しています。下のコードは _binary_clf_curve からいま最低限必要な動作のみ抜き出し、コメントをアレンジしたものです。

import numpy as np
from sklearn.utils.extmath import stable_cumsum

def my_binary_clf_curve(y_true, y_score, pos_label=None):
    y_true = np.array(y_true)
    y_score = np.array(y_score)

    # y_true および y_score を y_score の降順にソートする
    desc_score_indices = np.argsort(y_score, kind="mergesort")[::-1]
    y_score = y_score[desc_score_indices]
    y_true = y_true[desc_score_indices]
    print("y_score(ソート済み): ", y_score)
    print("y_true(ソート済み): ", y_true)

    # ほしいのはこの y_score のそれぞれを閾値としたときに、検出されるスパムメールと非スパムメールの個数のベクトル
    # ただ y_score はしばしば同値データを多く含む
    # 次の要素との差がゼロでないインデックスでのみスパムメールと非スパムメールの個数がほしいのでそのようなインデックスを取得する
    distinct_value_indices = np.where(np.diff(y_score))[0]
    print("次の要素との差がゼロでないインデックス: ", distinct_value_indices)
    threshold_idxs = np.r_[distinct_value_indices, y_true.size - 1] # 最後の fp、tp が取れないので最後のインデックスは補う

    # y_true を cumsum すればよい    
    tps = stable_cumsum(y_true)[threshold_idxs] # その閾値で検出されるスパムメールの個数
    fps = 1 + threshold_idxs - tps # その閾値で検出される非スパムメールの個数(インデックス+1 が検出されるデータ総数なので tps を引けばよい)
    return fps, tps, y_score[threshold_idxs]    

print("y_score: ", list_score)
print("y_true: ", list_label)
fp, tp, thresholds = my_binary_clf_curve(list_label, list_score)
print("fp: ", fp)
print("tp: ", tp)
print("thresholds: ", thresholds)
y_score:  [2, 1, 2, 4, 2, 1, 3, 5]
y_true:  [0, 0, 0, 1, 1, 0, 1, 1]
y_score(ソート済み):  [5 4 3 2 2 2 1 1]
y_true(ソート済み):  [1 1 1 1 0 0 0 0]
次の要素との差がゼロでないインデックス:  [0 1 2 5]
fp:  [0. 0. 0. 2. 4.]
tp:  [1. 2. 3. 4. 4.]
thresholds:  [5 4 3 2 1]

fp、tp は閾値を thresholds にしたがって変化させていったときに検出される非スパムメールスパムメールの個数です。呼び出し元の roc_curve ではこの fp と tp をそれぞれ最後のインデックスの値で割って割合(fpr、tpr)に直しています。drop_intermediate=True の場合は先に不要な座標の除去が行われます。この fpr、tpr が sklearn.metrics.auc に渡されます。

そしてこの auc から必要な動作を抜き出すと以下です。というか numpy.trapz に面積を求めさせているだけです。つまり、ROCカーブの下の面積は台形則で求められている(各座標を直線で結んだ折れ線の下の面積が求められている)ことがわかります。

def my_auc(x, y):
    area = np.trapz(y, x)
    return area

fpr = fp / fp[-1]
tpr = tp / tp[-1]
auc = my_auc(fpr, tpr)
print(auc)
0.9375

いかがでしたか? AUCの計算方法はとても簡単でしたね。やはり自分で適当に実装できそうです。
何がいいたいのかというと自分は適当にやって失敗したんですか、適当だからといってやるべきではないのは、スコア降順にソートした上でROCカーブの座標を1点ずつ拾っていくことだと思います。これをやるとROCカーブはスコア同値のデータの中でラベル1とラベル0がどう並んでいるかに依存します。「スコア同値ならばラベル昇順」というソートにすれば台形則から三角形を削った長方形則(というのか?)のAUCは出るといえば出ると思います(ただAUCの定義はやはり「閾値を下げていったときの fpr-tpr 曲線」の面積なので、定義に沿っていないかなり粗い下からの評価になってしまうと思います)。この方式だと三角形が削れているので scikit-learn より小さいAUCになります。じゃあこれに三角形を補えばいい気もしますが、カーブの座標からはどう三角形を補うべきかの情報が失われています(長方形の短冊を埋めるような三角形を補う、だと誤りです。同じスコアによって描かれた辺だけを三角形で埋める必要があります)。最初から素直に scikit-learn と同じ方法でやればいいと思います。
スパムメールを学習するモデルというのは普通は未知のメールに対してスパムメールか否か(1 or 0)を判定することを目指すものだと思います。なのに、ここでのモデルはスパムメールらしさのスコアというものを出力していて、そのスコアが AUC で評価されているのは不思議さがあると思います。まず、スパムメールか否かを判定するモデルでは、モデルはメールの色々な特徴からスパムらしさを説明しようとすると思います。モデルが最終的に実際に利用されるときは、色々な特徴から推測されるスパムらしさの合計値をどこかしらの閾値で切ってモデルの出力を "1 or 0" に丸めると思いますが、モデルの内部では "1" と判定したデータの中にも「スパムらしい特徴をとても強くもつなあ」「スパムらしい特徴をまあまあもつなあ」という違いはあると思います。AUC はその 1 or 0 に丸められる前のモデルのスパムらしさのスコア付けを評価するものです。なぜ 「(閾値をどこかしらに決めた)1 or 0 の最終判定における精度や感度」などではなく「(閾値を決める前の)スコアの AUC」をみたいのかというと、「そもそも実は最終判定ではなくスコアがほしい」というケースがあると思います。閾値は後から適当に調整するので、スコアをよこせということです。あるいは、元より最終判定のみがほしい場合でも、精度・感度は同程度のモデルが2つ出来上がってきたとして、2つとも検出漏れはないが誤検出はあるとして、スパムと誤判定された非スパムメールに対する内部スコアが「閾値は超えてしまったが少し超えてしまっただけで多くのスパムメールのスコアよりは小さい」というモデルと、「多くのスパムメールのスコアよりも大きくなってしまっている」というモデルだったら後者の方が判断の根拠が心配な気はします。この場合、後者の方がAUCが小さくなります。

論文読みメモ: GBrank(その1)

以下の論文を読みます。

Zhaohui Zheng, Keke Chen, Gordon Sun, Hongyuan Zha. A Regression Framework for Learning Ranking Functions Using Relative Relevance Judgments. In SIGIR, pages 287-294, 2007. https://www.cc.gatech.edu/~zha/papers/fp086-zheng.pdf
※ キャラクターは架空のものです。解釈誤りは筆者に帰属します。何もわかっていないのでおかしな点はご指摘ください。

f:id:cookie-box:20190101155733p:plain:w60

「相対的な関連度判定を用いた、ランキング関数学習のための回帰のフレームワーク」ですか。ランキング関数とは何でしょうか。

f:id:cookie-box:20190101160814p:plain:w60

アブストラクト冒頭に「商用検索エンジンにとって、効果的なランキング関数は不可欠だ」とあるね。ここでいうランキング関数とは、「ユーザが検索したフレーズを受け取って、提示可能な全てのウェブサイトにその検索フレーズとの関連度(relevance)を割り当てる」ような関数みたいだね。検索フレーズ-ウェブサイト対を受け取ってその関連の強さを返す関数といった方が入出力がはっきりするかな。検索エンジンの向こう側では何らかのランキング関数がユーザの検索フレーズと全てのウェブサイトとの関連度を計算して、その関連度が大きい順にウェブサイトが検索結果画面に並べられて提示されるわけだ。例えば「チーズタルト」とGoogle検索すると、1位にクックパッドのレシピ検索結果、2位にチーズタルト専門店のサイト、3位にまた別のチーズタルト専門店のサイトが表示されるけど、これはGoogle検索が使用しているランキング関数が「『チーズタルト』とレシピサイトの関連度 > 『チーズタルト』と専門店Xサイトの関連度 > 『チーズタルト』と専門店Yサイトの関連度」と関連度を算出したんだね。

f:id:cookie-box:20190101155733p:plain:w60

なるほど…私たちが普段検索している裏ではそんなことが…想像したことがありませんでした。しかし、そのようなランキング関数をどうやって学習するんでしょう。

f:id:cookie-box:20190101160814p:plain:w60

まず全ての「検索フレーズ-ウェブサイト対」は特徴ベクトルで表現されるみたいだね(論文3.1.1節)。その特徴ベクトルは検索フレーズのみに依存する特徴ウェブサイトのみに依存する特徴検索フレーズとウェブサイトの両方に依存する特徴を含むらしい。例えば以下だと挙げられているね。あくまで例だとは思うけど、でもこんな感じの特徴をつかっているんだろうね。

  • 検索フレーズが何単語からなるか
  • 検索フレーズが人名を含むか
  • そのウェブサイトに張られているリンクがいくつ存在するか
  • そのウェブサイト内にアンカーテキストか何バイトあるか
  • そのウェブサイトが何語のサイトか
  • そのウェブサイト内に検索フレーズ内の各単語が何回登場するか
  • そのウェブサイト内のアンカーテキストに検索フレーズ内の各単語が何回登場するか
ただこれらの特徴を用意しても、何を目指すべきかがないと学習はできない。この「検索フレーズ-ウェブサイト対」とこの「検索フレーズ-ウェブサイト対」ではどちらを上に表示するべきか、とかの情報がほしいよね。これについては、「人手で付けた優先度ラベル」を用いるやり方と、「ユーザの実績クリックデータ」を用いるやり方があるみたいだね。「人手で付けた優先度ラベル」はそのままの意味で、人手で「検索フレーズ-ウェブサイト対」に 0, 1, 2, 3, 4 のグレードを付けていくものらしい。ランキング関数がやるべきことは、よりグレードの高い「検索フレーズ-ウェブサイト対」により高い関連度を出力することだね。「ユーザの実績クリックデータ」の方は、実際の検索結果のログからクリック率(クリックされた回数 ÷ その検索結果が表示された回数)に有意差があると認められる「検索フレーズ-ウェブサイト対」のペア(対のペアってややこしいね)を抽出するらしい。ランキング関数はそれらのペアに対してクリック率が高い方により高い関連度を出力することを目指して学習する。

f:id:cookie-box:20190101155733p:plain:w60

人手でラベリングとは…人海戦術なんでしょうか。「そのウェブサイト内に検索フレーズ内の各単語が何回登場するか」などはわかりやすい特徴ですね。「チーズタルト」と検索した人が求める情報があるサイトには、「チーズタルト」と書いてあるでしょうし。サイト内のアンカーテキストというのも特徴として重視されているんですね。アンカーテキストがたくさんあるサイトは「これについて知りたい人はこちら」と道案内してくれる、優良なサイトということなんでしょうか。…なんと、数々のSEO対策記事にアンカーテキストの重要性が説かれていますね。このブログのアクセス数が伸びないのは、アンカーテキストが足りないということだったんですね…アンカーテキストを増やして検索順位を上昇させなければ…。

f:id:cookie-box:20190101160814p:plain:w60

やめて! このブログのアクセス数が伸びないのは単純に意味わかんないからだよ!

f:id:cookie-box:20190101155733p:plain:w60

この論文でやりたいことはわかってきた気がします。でも、例えばその人手でラベリングした 0, 1, 2, 3, 4 のグレードに合致したランキングになるように関連度を出力したいというのは、単純にその 0, 1, 2, 3, 4 を目的変数にして特徴ベクトルで回帰すれば済む話ではないでしょうか。特別に「ランキングの学習」なるものが必要なんでしょうか。

f:id:cookie-box:20190101160814p:plain:w60

それはどうかな。私も詳しくないしまだこの論文の中身を読んでないからわかんないけど、「0, 1, 2, 3, 4 を目的変数にして誤差2乗和が最小になる関連度を出力する」と、「0, 1, 2, 3, 4 のグレードと大小関係が合致するような関連度を出力する(グレード0からグレード4が等間隔に並ぶような関連度がほしいわけではない)」というのは少し違う話な気がするよ。それに、「チーズタルト」という検索フレーズの場合の「0, 1, 2, 3, 4」と「パンケーキ」という検索フレーズの場合の「0, 1, 2, 3, 4」は同じ座標系でなくてもいいし。ディープラーニングで座標変換ごと学習するなら一緒かもしれないけどね。

f:id:cookie-box:20190101155733p:plain:w60

確かに。…となると、評価指標がわからなくなってきました。どのようなランキング関数であればよいランキング関数といえるのでしょうか。出力する関連度とグレードとの誤差2乗和が小さい、というわけにはいきませんよね?

f:id:cookie-box:20190101160814p:plain:w60

論文の3.2節が「評価指標」だね。この論文では以下の3つの指標を用いるとあるよ。

  • Number of contradicting pairs: 優先順位付けを誤ってしまった(関連度の大小を正しい優先順位と逆にしてしまった)「検索フレーズ-ウェブサイト対」ペア数だね。
  • Precision at K%: まず全ての「検索フレーズ-ウェブサイト対」ペアに対してランキング関数が出力する関連度の差の絶対値を計算する。この絶対値が大きい方からK%のペアたちを取ってきて、優先順位付けの正解率をみる。それが Precision at K% みたいだね。これは、「ランキング関数が特に自信をもって関連度が違うと判定しているペアたち」に絞って判定の正確さをみるといった感じの指標だね(自信をもって、というのは誤解を招くかもしれない表現だけど)。
  • Discounted Cumulative Gain (DCG): これは真の適合度(人手でラベリングしたグレードと似ているけど、0, 1, 2, 3, 4 じゃなくて 0, 0, 3, 7, 10 みたいにもっとコントラストを付けたものかな。人間の実感に合った重みにしたものだと思う。グレードと対応しているとも限らないかもしれない)を何か決めて、各検索フレーズに対する検索結果内で真の適合度が大きいウェブサイトをより上に並べたかどうかをみる指標だ。具体的には、検索結果のウェブサイトの真の適合度を足していくんだけど、1位は一番大きい重みで、2位はもう少し小さい重みで、3位はさらに少し小さい重みで…というように、重みを変化させながら足す(具体的には、i 位のウェブサイトの真の適合度は  1/ \log_2 (i+1) の重みで足す)。これが DCG だ。ばっちり真の適合度の降順に並べていれば最大の DCG が達成される。「よいウェブサイトをより上の方に並べたか」ともいうべき指標だね。
この内最後の DCG は「絶対的優先度」が所与な場合にしか計測できない。つまり、人手で決めたグレードとか真の適合度とかがある場合はそれを利用して DCG が計算できるけど、ユーザの実績クリックデータから抽出したクリック率に有意差がある検索フレーズ-ウェブサイト対のペアしかないような場合は、ぺア間の「相対的優先度」があるだけで、「絶対的優先度」はわからないから、DCG は計算できないんだよね。

f:id:cookie-box:20190101155733p:plain:w60

3つも指標があるんですか? ランキングのよさの指標とは一本化できないものなんですね…。DCG が最も総合的な指標である感じもしますが、常には計測できないと。であればペア間の相対的優先度をどれだけ正しく当てたかくらいしか計測できませんが、しかし誤ったとしてもランキング関数が予測した関連度は僅差だったかもしれません。それも考慮するならば、Precision at K% がほしくなりますね。…うーん、やはり複数の指標があった方がよさそうです。でも何となくわかってきましたよ。おそらくランキング関数はこの「ペア間の相対的優先度をどれだけ正しく当てたか」を目的関数に学習するのでしょう? 論文タイトルにも「相対的な」とありましたし。ですから、相対的優先度が判明している全てのペアに対してペアの1つ目が優先されるなら +1、ペアの2つ目が優先されるなら -1 となるように2値分類すれば…あ、いえ、ほしいのはペアを入力して優先度を判定するモデルではありません。ペアの片割れだけを入力して、それに対して関連度スコアを出力してくれるモデルでなければ…いったいどうすれば…。

f:id:cookie-box:20190101160814p:plain:w60

少しテクニカルだからね。この論文では関連度を学習するのに勾配ブースティング木をつかう。1時点前の強学習器が相対優先度の判定を誤ったペアに対して、1時点前の関連度の判定を取り換えっこして、取り換えっこした判定を目的変数に弱学習器を学習する。

f:id:cookie-box:20190101155733p:plain:w60

なるほど。そういわれれば確かに強学習器の誤りは是正されていきそうですが…しかし、そのような学習をすることは結局何を目指しているのでしょうか。

(次回があれば)つづく

雑記: 決定木の話(途中)

別の記事で勾配ブースティング手法を勉強しようとしているのですが、そこでよく用いられる決定木について知らなかったので決定木の話を書きました。初学者なのでおかしい点がありましたらご指摘いただけますと幸いです。
参考文献

  1. はじめてのパターン認識 | 平井 有三 |本 | 通販 | Amazon
    • 11章で決定木を扱っています。
  2. 統計的学習の基礎 ―データマイニング・推論・予測― | Trevor Hastie, Robert Tibshirani, Jerome Friedman, 杉山 将, 井手 剛, 神嶌 敏弘, 栗田 多喜夫, 前田 英作, 井尻 善久, 岩田 具治, 金森 敬文, 兼村 厚範, 烏山 昌幸, 河原 吉伸, 木村 昭悟, 小西 嘉典, 酒井 智弥, 鈴木 大慈, 竹内 一郎, 玉木 徹, 出口 大輔, 冨岡 亮太, 波部 斉, 前田 新一, 持橋 大地, 山田 誠 |本 | 通販 | Amazon
    • 9章で決定木を扱っています。
  3. Induction of Decision Trees
    • ID3 開発者による ID3 及びその改善点についての論文で1986年のものですが、ID3 自体は1979年に発表されたものであるようです。C4.5 につながる連続値の取り扱いや情報量ゲイン比についても記述されています。
  4. A Brief History of Classification and Regression Trees
    • 拾い物のスライドですが、決定木の歴史について詳しそうです。
  5. A comparative study of decision tree ID3 and C4.5
    • C4.5 について記述がありそうだったのでメモしましたが、ちゃんと読んでいません。
  6. ID3 - Wikipedia
    • 例を借用したんですが、実は ID3 自体はこの例のような3値以上の分類に対応していないようです。

Appendix1. 自分が初めて決定木を学ぶ上で勝手に混乱したという読み飛ばしていい話
決定木について調べると、「決定木生成アルゴリズムの代表例として ID3、C4.5、CART などがある」と紹介されていることが多いと思います。それらが決定木生成アルゴリズムであることに間違いはないですが、一つのアルゴリズムというよりはむしろ必要な要素アルゴリズムを組み合わせて実行できるようにしたプログラム/システム(それもかなり古い)です(のはず)。なので、それらについて調べようとしても明確な数学的定義というよりは当時こんなものが発表されたという断片的な情報が多いです(CART については当時原書籍が出たことしか自分はわからなかったので、そもそも各要素が一義的にこれと定められていたのかも自分はよくわかっていません)。また、ID3、C4.5、CART の要素アルゴリズムの組み合わせは必然性が強いというものではないはずです。なので、ID3、C4.5、CART などが決定木を学ぶ軸になるというものではないです。自分で学ぼうとするとき、それらを踏まえないと混乱すると思います。
そもそも ID3(1979年)と C4.5(1993年)は Quinlan さんが開発した分類木生成プログラムであり、CART1984年)は Breiman さんらがまとめた回帰/分類木生成システムです(のはず)(ただし CART といったとき、現在まで改修され続けている CART の血を継ぐ各種パッケージのこととか、もっと一般的な回帰/分類木を生成する方法論を指していることも多いと思います)。決定木を生成する典型的な方法は「分岐に利用する特徴と境界をどのように選ぶか」「どこまで分岐させるか」「終端にどんな予測値を割り当てるか」の3つの要素をどう実現するか決めることで(参考文献1. 178ページ)、ID3 も C4.5 も CART もこの3つの要素を何か実装しています(かつ、その実装により対応できる入出力の型がカテゴリ値のみか連続値も可能か決まっています)。ただそれはプログラムの実装がそうなっているというだけです。例えば ID3 は連続入力値を扱えないとされますが、いま私たちが再実装するならば閾値探索を入れる拡張は容易にできるはずです。なんというか、ID3 の特徴選択基準(情報量ゲイン)が絶望的に連続値に対応できないというわけではないです。
いいたいことは、ID3、C4.5、CART は厳密にこんな要素の組み合わせなのだと明らかにしようとするのは歴史の勉強でしかなくて(もちろんその組み合わせで概してどのような学習の傾向がみられたという情報は重要ですが)、大事なのは決定木生成アルゴリズムがどのような要素アルゴリズムからなり(=上述の緑太字)、それぞれの要素アルゴリズムにどのようなバリエーションがあり、それぞれのバリエーションがどのような意味をもつのかだと思います。無論、上の3つの要素の組み合わせにあてはまらない木もあります(階層的混合エキスパートなど)。
Appendix2. ID3、C4.5、CART のまとめ
Appendix1. にかいた通り ID3、C4.5、CART が厳密にどうなっているかを追いかけることは本質的ではないと思いますが、どのようなアイデアが含まれているか自体は重要と思うので自分の認識をまとめておきます。正確ではないです。
ID3C4.5CART
入力ベクトルの型カテゴリ値(欠損不可)カテゴリ値 or 連続値(欠損可)
目的変数の型2値カテゴリ値(3値以上の場合は各カテゴリ値かどうかで別々に木をつくる)連続値
分岐の仕方のよさの指標分岐後のノードたちのエントロピーのデータ数重み付き平均が小さいような分岐ほどよい分岐とする。
⇔ 分岐元ノードのエントロピーに比べて分岐後のノードたちのエントロピーの重み付き平均がどれだけ減ったかをその分岐の情報量ゲインと定義し、情報量ゲインが大きい分岐ほどよい分岐とする。
ID3 と似ているが、情報量ゲインは分岐が多いほど大きくなりがちなので、過学習を防ぐ観点から、情報量ゲインを分岐のエントロピー(たくさん枝分かれさせると大きい)で割った情報量ゲイン比に修正し、情報量ゲイン比が大きい分岐ほどよい分岐とする。 典型的には分岐後のノードたちのジニ指数のデータ数重み付き平均が小さいような分岐ほどよい分岐とする。ただノードにおける目的変数の不純度を測る指標なら何でもよくて、エントロピーなどでもよいはず。 分岐後の誤差2乗和が最小となる分岐ほどよい分岐とする(他の誤差の指標も適用できそうだが)。
その指標にしたがって分岐させる方法各カテゴリ特徴について、カテゴリ値が取りうる値だけ分岐させたときの分岐の情報量ゲインを計算し、最も情報量ゲインの大きい特徴で分岐させる。 特徴がカテゴリ値の場合は、C4.5はID3同様にカテゴリ値がとりうる値だけ分岐させたときの分岐のよさの指標を計算する? CARTについては2分岐しかしないという情報を散見するので、3値以上とりうるカテゴリ値はone-hotにするんだと思う。パッケージによると思う。
特徴が連続値の場合は閾値によって2つの枝に分岐させるすべてのパターンについて分岐のよさの指標を計算する。
何にせよ、最も分岐のよさの指標がよい特徴および分岐パターン/閾値で分岐させる。
一旦どこまで分岐させるかのルールノードに同一ラベルのデータしか存在しなくなったらそこで止める。または、分岐させることのできる特徴がなくなってしまったらそこで止める。または空のノードができたらそこで止める(とりうるカテゴリ値の枝を全てつくるらしいので空ノードがありうる)。参考文献2. の記述では、終端ノードに振り分けられるデータ数がある値まで小さくなるまで分割するというようにある(351ページ)。
一旦どこまで分岐させた後に余分な枝を刈り取るルール刈り取らない。刈り取る。後で追記。「判定の誤りの指標」に「終端の数」を適当なバランスで足し合わせた「修正誤り指標(木のコスト)」をつくって、修正誤り指標の削減に貢献していない枝を刈り取る。
終端ノードに割り当てる予測値通常は終端ノードに振り分けられたデータの目的変数の平均値。後で追記。

決定木とは

決定木はあなたにまず最初の質問を出してきます。あなたがその質問に答えるとその答えによって別の異なる質問を出してきます。その質問にも答えるとまた別の異なる質問を出してきます。それを何度か繰り返したのち、あなたに何か判定結果を下してきます。そんな木です。要するにアキネイターです。

やりたいこと

既にいるアキネイターから判定結果を得るには質問に答えていけばいいです。いまやりたいのは新しいアキネイターをつくることです。言い換えると、
  • 目的は個々の入力データに対して何らかの判定をすることです(※1)
  • それを達成する方針として、入力データに、入力データが答えらえる質問を繰り返すことにします。
  • やるべきことはこの「どんな質問たちを、どんな順番で出して、どんな判定を下すか」のマニュアルを完全に決めることです。そういうマニュアルを学習データを利用してつくります。
完成したマニュアルはきっと、最初の質問が書いてあって、そこから回答によって枝分かれした先のそれぞれにまた質問が書いてあって、また枝分かれして…終端には判定が書いてあるものになっていると思います。まるで木です。
※1. 判定するとは例えば、服屋に来店する個々のお客さんがどの商品を欲しがるかとかです。服屋に来店するお客さんという入力データには「男性か女性か」「子どもか大人か」「背が高いか低いか」のような質問ができそうです。「江戸時代のご先祖が武士だったか」のような、来店するお客さんのデータで答えられない質問は駄目です。また、「どの商品が欲しいか」自体を訊くのもいまは駄目です。というかそれが訊けるなら訊いた方がいいですが、新しく来店したお客さんはまだこの店のラインナップをよく知らないので、こちらで何がお気に召しそうか予測して売り場にご案内したいという状況とでも考えてください。

やりかた


しかしそのようなマニュアル=木をつくるのは大変そうです。入力データに訊くべき質問は色々ありえるし、質問の順番も決めなければなりません。学習データに対してなるべく正しい判定をする木を探せばよさそうですが、極端な話、何も考えずに木をどこまでも枝分かれさせれば個々の学習データを完全に正しく判定する木というのができてしまいそうです。しかし、そんな木はきっと学習データの思わぬノイズまでも学習してしまって未知データへの性能が悪そうです。なので、「ノイズらしきデータからはなるべく学習しない」を満たすような木のつくり方をしたい気がします。
基本的な決定木生成アルゴリズムでは、質問は以下のようにすることに決めています。
  • 1回の質問では「1つの特徴の値がどんなか」のみを訊く(※2)
そして、最初に学習データをなるべくよく振り分ける質問をして、それで振り分けたそれぞれのデータに対してまたなるべくよく振り分ける質問をして…を繰り返して木を生成します。後で追記。


※2. つまり「2つの数値的特徴の積がどんなか」のような質問は訊かないことにしているわけです。アルゴリズムの中でそのような質問も考慮してみてもいいかもしれませんが、それを常に調べるのはとても時間がかかりそうですし、個別の特徴への質問を組み合わせて繰り返せば近似的に「2つの特徴の積がどんなか」っぽい振り分けは達成できそうです。しかし、ある2つの特徴の積を考慮した方がよいと信じられる状況であれば、それを新たな特徴として追加しておくのがよいと思います。

誤差2乗和による回帰木の生成

後で追記。

情報量ゲインによる分類木の生成

後で追記。下のスクリプトを参照。

ジニ指数による分類木の生成

後で追記。

計算用スクリプトアルゴリズムの実装ではありません)

情報量ゲインによる分類木の生成です。ウィキペディアのID3の例を借りていますが、実はID3自体はこの例のような3値以上の分類に対応していません。目的変数が3値以上でもエントロピーを計算することはできますが、「それでエントロピーが小さい方がよりよい分岐の仕方といえるのか」というところがあやしくなってくるんだと思います。しかし下の例ではそういう細かい点は無視することにしてアルゴリズムを実行しています。

import numpy as np
import pandas as pd
from collections import OrderedDict

# 準備.pandas.Series を渡すとユニーク要素ごとの出現確率のベクトルを返す関数
def get_probabilities(series):
    probabilities = series.value_counts() / series.shape
    return probabilities

# 準備.離散確率分布を渡すと平均情報量を返す関数(ウィキペディアに合わせて対数の底=3)
def get_entropy(probabilities):
    return np.sum(probabilities * np.log2(1.0 / probabilities) / np.log2(3.0))

# 準備.データとラベル名と特徴名とその特徴のあるカテゴリ値を渡すと
# 「各特徴が各カテゴリ値か否か」で2ノードに分岐させたときの
# 各ノードの平均情報量を各ノードに振り分けられるデータ数で重み付き平均した値を返す関数
def get_expected_entropy(df, label_name, feature_name, feature):
    df0 = df[df[feature_name]!=feature]
    w0 = df0.shape[0] / df.shape[0]
    df1 = df[df[feature_name]==feature]
    w1 = df1.shape[0] / df.shape[0]
    return w0 * entropy(get_probabilities(df0[label_name])) \
           + w1 * entropy(get_probabilities(df1[label_name]))

# ------------------------- やりたいことはここから -------------------------
# いまやりたいことは、以下の「食性」「発生」「体温」の特徴から「分類」を予測すること
od = OrderedDict()
od['名前'] = ['ペンギン', 'ライオン',   'ウシ', 'トカゲ', 'ブンチョウ']
od['食性'] = [    '肉食',     '肉食',   '草食',   '肉食',       '草食'] 
od['発生'] = [    '卵生',     '胎生',   '胎生',   '卵生',       '卵生'] 
od['体温'] = [    '恒温',     '恒温',   '恒温',   '変温',       '恒温'] 
od['分類'] = [    '鳥類',   '哺乳類', '哺乳類', '爬虫類',       '鳥類'] 
df = pd.DataFrame(od)
print(df)

# まずすべてのデータを根ノードに所属させる
# 根ノードでの平均情報量をみる
ent = get_entropy(get_probabilities(df['分類']))
print(ent) # 0.9602297178607613 ... 平均情報量が大きいほどこのノードでの混じり気が大きい

# 混じり気が小さくなるように分岐させたいので、「各特徴が各カテゴリ値か否か」で2ノードに分岐させたときの
# 各ノードの平均情報量を各ノードに振り分けられるデータ数で重み付き平均した値が、
# いまの混じり気よりもどれだけ減るかをみる
print(ent - get_expected_entropy(df, '分類', '食性', '肉食')) # 食性が肉食か否か --> 0.1078578164321784
print(ent - get_expected_entropy(df, '分類', '発生', '卵生')) # 発生が卵生か否か --> 0.6126016192893443
print(ent - get_expected_entropy(df, '分類', '体温', '恒温')) # 体温が恒温か否か --> 0.45548591500359525

# 「発生が卵生か否か」が最も混じり気を減らすのでそれで分岐させることに決める
df_0 = df[df['発生']!='卵生']
df_1 = df[df['発生']=='卵生']
ent_0 = get_entropy(get_probabilities(df_0['分類']))
ent_1 = get_entropy(get_probabilities(df_1['分類']))
print(ent_0) # 0.0 --> こちらのノードではすでに混じり気なし(発生が胎生なら分類は哺乳類でしかない)
print(ent_1) # 0.5793801642856949 --> このノードではまだ混じり気あり

# 混じり気を小さくしたいので、また分岐を試みる
print(ent_1 - get_expected_entropy(df_1, '分類', '食性', '肉食')) # 食性が肉食か否か --> 0.15876032857138994
print(ent_1 - get_expected_entropy(df_1, '分類', '体温', '恒温')) # 体温が恒温か否か --> 0.5793801642856949

# 「体温が恒温か否か」が最も混じり気を減らす(というか完全に混じり気を減らす)のでそれで分岐させることに決める
df_1_0 = df_1[df_1['体温']!='恒温']
df_1_1 = df_1[df_1['体温']=='恒温']
ent_1_0 = get_entropy(get_probabilities(df_1_0['分類']))
ent_1_1 = get_entropy(get_probabilities(df_1_1['分類']))
print(ent_1_0) # 0.0 --> 混じり気なし(卵生かつ変温なら爬虫類)
print(ent_1_1) # 0.0 --> 混じり気なし(卵生かつ恒温なら鳥類)

# よって、やりたいことは、以下の木にしたがって分類することで達成される
# 「発生が卵生か」 -- No --> 哺乳類 
#                  -- Yes --> 「体温が恒温か」-- No --> 爬虫類
#                                             -- Yes --> 鳥類

NeurIPS2018読みメモ: CatBoost: unbiased boosting with categorical features(その-1: そもそも gradient boosting がわからない)

以下の論文を読みます。

Liudmila Prokhorenkova, Gleb Gusev, Aleksandr Vorobev, Anna Veronika Dorogush, Andrey Gulin. CatBoost: unbiased boosting with categorical features. In Proceedings of NIPS 2018. https://papers.nips.cc/paper/7898-catboost-unbiased-boosting-with-categorical-features
※ キャラクターは架空のものです。解釈誤りは筆者に帰属します。何もわかっていないのでおかしな点はご指摘ください。

f:id:cookie-box:20190101155733p:plain:w60

CatBoost…直訳すると「猫が励ます/後押しする」…なるほど、ツイッターなどをみると猫に励まされている方が大勢いますものね。

f:id:cookie-box:20190101160814p:plain:w60

違うよ?

f:id:cookie-box:20190101155733p:plain:w60

はい。この論文は CatBoost という新しい「勾配ブースティング(gradient boosting)」のツールについて解説するものであるようです。こ、勾配ブースティングとは? 勾配が励ますということですか? 猫ちゃんが励ましてくれるのではなく?

f:id:cookie-box:20190101160814p:plain:w60

いや、機械学習の文脈なんだから、「勾配ブースティング」というのは、モデルを理想のモデルに近づけるためにどんどん更新していくときに損失関数が小さくなる方向(勾配方向)に何か効果的に更新できる手法ということじゃないのかな…あ、ウィキペディアに「勾配ブースティング」という記事があるよ。

これを読むと、勾配ブースティングは「予測モデルを複数個組み合わせて(アンサンブルさせて)最終的な予測を下す」ような学習をする場合の学習手法の一つらしい。このときの個々のモデルを「弱学習器」、組み合わせた全体のモデルを「強学習器」というみたいだね。

f:id:cookie-box:20190101155733p:plain:w60

「弱学習器」を組み合わせて「強学習器」にする、ですか…でも、適当にその辺の弱い人々を100人集めてチームを組ませたところで藤井聡太七段に将棋で勝つことはできませんよね?

f:id:cookie-box:20190101160814p:plain:w60

弱学習器に「訓練されていない」って意味はないからね!? ちゃんと訓練して!? というか弱学習器たちをどう訓練するべきかって手法が「勾配ブースティング」だからね。それでその勾配ブースティングでは、以下のように弱学習器 h_m(x) \; (1 \leqq m \leqq M) を足し上げて強学習器 F_M(x) にするらしい。各 m ステップでは1時点前の強学習器 F_{m-1}(x) の残差  y - F_{m-1}(x) を埋めるような弱学習器 h_m(x) を得ようとするイメージだと紹介されているね。

  •  \hat{y} = F_1 (x) \equiv h_1(x) で与えられる \hat{y} がよい予測値になるような h_1(x) を得る。
  •  \hat{y} = F_2 (x) \equiv F_1(x) + h_2(x) で与えられる \hat{y} がよい予測値になるような h_2(x) を得る。
  •  \hat{y} = F_3 (x) \equiv F_2(x) + h_3(x) で与えられる \hat{y} がよい予測値になるような h_3(x) を得る。
  • \cdots
  •  \hat{y} = F_M (x) \equiv F_{M-1}(x) + h_M(x) で与えられる \hat{y} がよい予測値になるような h_M(x) を得る。
もう少しきちんと説明すると、いま、考えている損失関数 L の平均を最小化するような関数 \hat{F} がほしい。
\hat{F} = \underset{F}{\rm arg \, min} \, E_{x,y} \bigl[L \bigl( y, F(x) \bigr) \bigr]
この \hat{F} をどう手に入れようかということになるけど、ここではある関数クラス \mathcal{H} に属する関数 h_m \in \mathcal{H}M 個もってきてその重み付き和(+定数)として得ることにする。
 \displaystyle \hat{F}(x) = \sum_{m=1}^{M} \gamma_m h_m(x) + {\rm const.}
それでこの \{h_m\} を得るのに、ここでは以下のように訓練データ上での損失を最小化することにしたい。
\displaystyle \begin{split} F_0(x) &= \underset{\gamma}{\rm arg \, min} \sum_{n=1}^{N} L(y_n, \gamma) \\ F_m(x) &= F_{m-1}(x) + \underset{h_m \in \mathcal{H}}{\rm arg \, min} \left[ \sum_{n=1}^{N} L\bigl( y_n, F_{m-1}(x_n) + h_m(x_n) \bigr) \right] \end{split}(☆)
けど、これで任意の損失関数 L について最適な \{h_m\} を得るのは一般に不可能なので、F_m(x) の更新式として上式の代わりに以下の式を用いることにする。これが勾配ブースティングだって。
 \displaystyle \begin{split} F_m(x) &= F_{m-1}(x) - \gamma_m \sum_{n=1}^N \nabla_{F_m-1} L \bigl(y_n, F_{m-1}(x_n) \bigr) \\ \gamma_m &= \underset{\gamma}{\rm arg \, min} \sum_{n=1}^{N} L \bigl(y_n, F_{m-1}(x_n) - \gamma \nabla_{F_{m-1}} L \bigl(y_n, F_{m-1}(x_n) \bigr) \end{split}
…思うに、本当は各ステップで最も上手く残差を埋める h_m を得たいけど(☆)、たぶん弱学習器がとりうる関数空間 \mathcal{H} 内でベストな h_m を見つけるのって難しい。そもそも ☆ は h_m を陽に示してないし。だから代替案として、「1時点前の強学習器 F_{m-1} での L の勾配方向 \nabla_{F_m-1} L \bigl(y_n, F_{m-1}(x_n) \bigr) 」を正解ラベルとして、\mathcal{H} 内でこれを最もよく当てる  h_m を見つけることにする。これだって厳密な最適解を見つけるのは難しそうだけど、いま「  \mathcal{H} から最もよく正解ラベルを当てる弱学習器を選ぶ」枠組みというのはあるはずだから(学習できる弱学習器を選んだはずだから)、それにあてはめればいい。問題は選んだ弱学習器  h_m を強学習器にどれだけ反映するかの係数 \gamma_m だけど、これは「1次元だから適宜損失が最小になるように探索して」って感じだね。そしたら強学習器を  F_m(x) = F_{m-1}(x) - \gamma_m h_m と更新すればいい。

f:id:cookie-box:20190101155733p:plain:w60

え、えっと? 要するに、2人目以降は前の人までの判断のミスを埋めることを目指して学習するんですね? そうですね…例えばスパムメール分類器を学習するとして、5つの訓練データの正解ラベルを (1,1,0,0,0) とでもしましょう。1人目が訓練データの特徴からこの正解ラベルを目指して学習した結果が (1,0,0,0,1) であったとします。この場合、2人目は1人目のミスを埋めるように学習しなければならないので、2人目にとっての正解ラベルは (0,1,0,0,−1) ですね、おそらく(符号はなんとでもなるのでどちら向きでもいいですが)。さらに3人目は1人目と2人目の判断を最適に合体させたような強学習器のミスを埋めることを目指すことになります…でもこれ、ミスを埋めていくといわれると確かに強学習器がどんどん改善されていきそうですが、1時点前までの強学習器の損失っぷりが「弱学習器で埋められるようなもの」でなければなりませんよね…? 例えば個々の弱学習器が線形回帰モデルで損失関数が2乗誤差の和だったら、1人目ができる限り頑張って学習した時点で残差は平均ゼロのノイズであるはずです。平均ゼロのノイズに線形回帰モデルをフィッティングしても傾きもバイアスもゼロですから、この残差を2人目が埋め合わせようというのはできないのでは…?

f:id:cookie-box:20190101160814p:plain:w60

1人目と2人目が参照する訓練データが全く同じなら損失が2乗誤差の線形回帰モデルではそうなりそうだね。でもロジスティック回帰による線形分類なら違いそうかな。ロジスティック回帰結果の真のラベルとの誤差を再びロジスティック回帰することはできそうだから。もっとも、こういう組み合わせて強くする学習に向いてるモデルの最たる例が「決定木」らしい。ウィキペディアの「勾配ブースティング」にも「勾配ブースティングは典型的には固定サイズの決定木を弱学習器として使用される」とあるね。

f:id:cookie-box:20190101155733p:plain:w60

決定木? 何ですかそれは? ご神木から得た神託を予測値とするのですか?

f:id:cookie-box:20190101160814p:plain:w60

違うかな。ウィキペディアに Decision tree learning という項目があるね。

上の記事の最初の図に決定木の例があるよ。"is sex male?" を根ノードにした木構造をしているね。これは「タイタニックの乗客が生存したか否か」の2値分類だ。各乗客の特徴として「性別」「年齢」「ともに乗船した配偶者または兄弟の人数」が与えられているらしい(この木で使用されていないだけで他の特徴もあるかもしれないけど)。各乗客は根ノードの質問から出発して、各ノードの質問の答えに応じていずれかの終端ノードに振り分けられる。緑色の終端ノードに振り分けられたら「生存」、赤色の終端ノードに振り分けられたら「死亡」と判定する。例えば女性の乗客なら真っ先に一番右の終端ノードに振り分けられて「生存」と判定されるけど、実際女性の乗客の生存率は 73% みたいだね。女性ならば「生存」といっておけば 73% は当たる。27% は外してしまうけど。というか図のたった3つの内部ノードだけで全乗客に対して 80% 以上の正解率が出るんだね( 0.36 * 0.73 + 0.02 * 0.89 + 0.02 * (1 - 0.05) + 0.61 * (1 - 0.17) = 0.8059 )。タイタニックの乗客で生存したのは雑にいうと「女性か、または9歳以下でかつともに乗船した兄弟が2人以下の男児」だったらしい。

f:id:cookie-box:20190101155733p:plain:w60

同じ9歳以下の男児でも兄弟の人数でこうも生存率が変わりますか…3人以上の兄弟と乗船した9歳男児は下の兄弟の避難を優先させたのか、あるいは三等船室の乗客であった確率が高かったのでしょうか…。しかし、「タイタニック」はもう21年前の映画なんですね。今時の高校生はタイタニックと聞いて「あー沈没した船ね」となるのでしょうか。

f:id:cookie-box:20190101160814p:plain:w60

…それで、この図のような決定木をつくる方法だけど、まず根ノードでは最も説明力の高い特徴の閾値によって分岐させて、そこから先は分岐させたデータの中で最も説明力の高い特徴の閾値によって分岐させて…といった手続きによるみたいだね。そうすればなるべくバランスのよい木はつくれるけど、データを決定付けるよい閾値を拾いきれるとは限らない。だから、いくつかの異なる木を組み合わせたくなってくる。組み合わせる方法の1つは訓練データの1部ずつをつかっていくつも木をつくって多数決させるというもので、この手法は「バギング」といってランダムフォレストもこの一種だ。組み合わせる別の方法は、バギングのように並列に木を組み合わせるのではなくて、直列に組み合わせる。それぞれの木はそれまでの木のミスを補うように学習する。これが「ブースティング」で、その具体的な強力なアルゴリズムの一つが「勾配ブースティング」だね。

f:id:cookie-box:20190101155733p:plain:w60

なるほど。先ほどの「勾配ブースティング」の手続きを、決定木を弱学習器として実行するわけですか。決定木であったとしても、各弱学習器 h_m はそれまでの強学習器の残差を正解ラベルとして学習し、得られた h_m を損失が最小になるような係数 \gamma_m で強学習器に足し合わせるというのは変わりませんね。

f:id:cookie-box:20190101160814p:plain:w60

なんだけど、決定木の場合はもう少し改良できるってかいてある。決定木って特徴空間を直方体の領域に分割するよね。この直方体ごとに \gamma_m を変えても構わない。というか、h_m の領域分割さえわかれば、h_m が各領域に割り当てた値に興味もなくて、各領域にどれだけの定数を足し合わせれば強学習器の損失が最小になるかを探索するだけだね。

f:id:cookie-box:20190101155733p:plain:w60

…確かに全ての領域について同じ係数で足してやらなければならないということもありませんね。新しい弱学習器は、特徴空間内の一部の領域の判定精度を向上させるが、その領域の外はむしろそれまでの強学習器のままでよく余計に何か足さないでほしい、ということもあるでしょうし。…しかし、こうウィキペディアだけ読んでいると雲をつかむようです。実際に何か決定木を学習してみることはできないでしょうか。

f:id:cookie-box:20190101160814p:plain:w60

決定木というかむしろ勾配ブースティング決定木のパッケージはいくつもあるから、何か動かしてみればいいんじゃないかな。

f:id:cookie-box:20190101155733p:plain:w60

そうですね、では xgboost というパッケージのR版をつかってみましょう(スクリプトは記事の下部)。データは…タイタニックのデータを拾ってくるのも面倒なので iris でいいです。木の数は指定できるのですね。しかし木を2つに指定しても6つできてしまうようですが…ああ、3クラス分類だから弱学習器1つにつき木が3つなのですね。面倒なので2クラスに絞りましょう。しかし、2つ目の弱学習器によって1つ目の弱学習器単体よりも分類精度がよくなる特徴と木の深さの組み合わせがなかなかないですね…特徴をがく片の幅と長さだけに絞り、木の深さを4にすると申し訳程度に訓練誤差が減少しますね。

[1]	train-error:0.260000 
[2]	train-error:0.240000 
   pred
y    0  1
  0 34 16
  1  8 42
        Feature       Gain     Cover Frequency
1: Sepal.Length 0.94330762 0.7858799 0.6666667
2:  Sepal.Width 0.05669238 0.2141201 0.3333333

f:id:cookie-box:20190201004758p:plain:w640
しかし木まで図示してくれるのですか。便利なものですね。

f:id:cookie-box:20190101160814p:plain:w60

こうなると特徴空間の分割もみたいよね。

f:id:cookie-box:20190101155733p:plain:w60

特徴空間の分割ですか…これでどうでしょう。白丸が "versicolor"、黒丸が "virginica" のデータで、破線が各学習器の分割境界、青が濃いほど "versicolor" と判定、赤が濃いほど "virginica" と判定という意味です。赤と青の濃さは上のグラフの Value に比例しています。一番左の図が1つ目の弱学習器、真ん中の図が2つ目の弱学習器で、一番右の図は結果の強学習器です。2つの弱学習器の図を単純に重ね合わせたら赤と青が相殺しないのがいけてないですね…。

f:id:cookie-box:20190201005400p:plain

f:id:cookie-box:20190101160814p:plain:w60

それどうやって描いたの??

f:id:cookie-box:20190101155733p:plain:w60

上のグラフを見て自力で長方形を描きましたね…。

f:id:cookie-box:20190101160814p:plain:w60

そっか…。

(次回があれば)つづく

スクリプト
library(xgboost)
max.depth <- 4 # 木の最大深さ
nrounds <- 2 # 木の数

### iris データの準備と特徴カラムと損失関数の定義
df <- data.frame(iris)
# いきなり3クラス分類だとわかりづらいので setosa を除外して2クラス分類にする
target_labels <- c("versicolor", "virginica") # "setosa"
num_class <- 1 # 3 # 3クラス以上でない場合は1を指定する
objective <- "binary:logistic" # "multi:softmax" 
df <- df[df$Species %in% target_labels, ]
df$Species <- factor(df$Species, levels=target_labels) # クラスを絞ったので factor 型を修正
# 特徴も2つのみに絞る
feats <- c("Sepal.Length", "Sepal.Width") # "Petal.Length", "Petal.Width"


# 説明変数と目標ラベル
x <- as.matrix(df[,feats])
y <- as.integer(df[,"Species"]) - 1

# 訓練する
model <- xgboost(data=x, label=y, max.depth=max.depth, nrounds=nrounds, 
                 objective=objective, num_class=num_class, verbose=1)

# 訓練したモデルの判定結果(訓練データ上)
pred <- predict(model, x)
if (length(target_labels) == 2) { # 2クラス分類にした場合は予測結果が 0~1 なので閾値で 0 or 1 に
    pred <- ifelse(pred < 0.5, 0, 1)
}
print(table(y, pred)) # 混同行列

# 各特徴の重要度を取得する
imp <- xgb.importance(feats, model)
print(imp) # 重要度

# モデルを図示
library(DiagrammeR)
print(xgb.plot.tree(feature_names=feats, model=model))
print(xgb.dump(model))

NeurIPS2018読みメモ: Deep State Space Models for Time Series Forecasting(その1)

以下の論文を読みます。

Syama Sundar Rangapuram, Matthias Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, Tim Januschowski. Deep State Space Models for Time Series Forecasting. In Proceedings of NIPS 2018. https://papers.nips.cc/paper/8004-deep-state-space-models-for-time-series-forecasting

※ キャラクターの原作とは関係ありません。解釈誤りは筆者に帰属します。問題点がありましたらご指摘ください。

f:id:cookie-box:20180305231302p:plain:w60

この論文のタイトルは「時系列予測のためのディープ状態空間モデル」だって。ちなみに NeurIPS2018 にはこの論文とは別に「条件なし単語生成のためのディープ状態空間モデル( Deep State Space Models for Unconditional Word Generation )」っていうタイトルめっちゃかぶってる論文もあるから注意な。

f:id:cookie-box:20180305232608p:plain:w60

時系列予測? そんなのどんな時系列かによるでしょう? 検証に利用しているデータは何ですか?

f:id:cookie-box:20180305231302p:plain:w60

なんで俺が説教されるの!? えっと、この論文では定性評価のための検証と定量評価のための検証とやってるみたいだけど、定性評価パートでは人工データのパラメータを復元できるかを検証してるみたいだから、現実の時系列データをつかってるのは定量評価パートだな。定量評価パートでつかってるのは、

  • "electricity": これは370人の顧客の時間ごとの電力消費量?
  • "traffic": サンフランシスコの高速道路の963の車線の時間ごとの占有率? 混み具合ってこと?
どっちも1日の中での「季節性」があるってかいてある。夜の方が電気がよく消費されるとか、通勤時間は道路が混んでるとか、そんな感じかな。public なデータみたいだから、ちょっとみてみるか。どこで手に入るんだ? まず以下の論文への参照があるな。なんかこの論文もこの論文で内容気になるけど、この論文が "electricity" と "traffic" データを検証につかってるみたい。で何々… "The data sets electricity and traffic are obtained from the UCI repository" ってあるから UCI リポジトリってところから取ってくるのか。じゃあ最初の論文でそうかいてくれよ…。とりあえず "UCI repository" で検索すれば出てきそうだな。
  • https://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014
    • "electricity" の方はたぶんこれだな…時刻のカラムを除くとカラム数が 370 ある。370人の顧客の電力消費量って370人の顧客の電力消費量の合計じゃなくて個々の顧客の電力消費量って意味だったのか…。数値は15分ごとの電力消費量(kW)で、時刻はポルトガル時間で毎年3月と10月にそれぞれ23時間しかない日と25時間ある日があるって。サマータイムだな。
  • https://archive.ics.uci.edu/ml/datasets/PEMS-SF
    • "traffic" の方はこれか。これは 10 分ごとの 963 個のセンサーの測定値みたいだ。PEMS_train ファイルは行ごとに 963 × 144 の行列になってるから、これが1日分のデータ(963個のセンサーの10分ごとの測定値=963×6×24=963×144)で、PEMS_trainlabels ファイルは PEMS_train ファイルの各行が「何曜日のデータだったか」を示してるな(1~7の各数字が月~日曜日を意味する)。曜日を識別する分類問題を意識したデータみたいだ。データ内には深夜にセンサーが止まって数値が欠損している日も2日含まれてるらしい…。
これらのデータの一部をプロットしてみるとこうかな。顧客370人分は多すぎるから5人分だけで、センサー963個は無理だから3個だけで、期間もほんの一部に絞ってあるけど。

f:id:cookie-box:20190121200857p:plain
f:id:cookie-box:20190121210707p:plain
f:id:cookie-box:20180305231302p:plain:w60

"electricity" は人が寝静まってる深夜の時間帯が落ち込んでる感じなのかな? みんな寝てるからパソコンしないとか? "traffic" はよくわかんないけど、やっぱり深夜のが道路が空いてる? 日曜日と比べて水曜日は朝方に道路が混んでるな。でも夜は若干日曜のが混んでる? って考えるとなんか特徴はありそうだけど、全部の曜日のパターンを見抜くのは難しそうだな…。ただこの論文では曜日の分類はしてなくてあくまで数値の予測をしようとしてるんだと思うけど。

f:id:cookie-box:20180305232608p:plain:w60

どちらも10~15分毎に観測されるデータで、異なる顧客/車線から同種のデータが大量に得られるんですね。

f:id:cookie-box:20180305231302p:plain:w60

アブストラクトによると、時系列予測のために状態空間モデルと深層学習を組み合わせたみたいなんだけど、時系列ごとの線形状態空間モデルを RNN で結合してパラメタライズすることで??、状態空間モデルと深層学習のいいとこどりをしてるって…状態空間モデルから受け継いだいいところは、解釈性と data efficiency?

f:id:cookie-box:20180305232608p:plain:w60

そのまま読めばデータを効率よくつかえるということなんでしょうが…カルマンフィルタのゲインのようなイメージでしょうか。カルマンゲインによって新しく観測された情報を必要なだけ取り入れるような。

f:id:cookie-box:20180305231302p:plain:w60

まあ詳しくは後で出てくるだろ。深層学習の方のいいところは複雑なパターンも学習できると。それで提案手法は、学習データが少ししかないときにデータの構造を想定するのにも、学習データが大量にあるときに複雑なパターンを学習するのにも効果的らしい…。Conclusions をみてもかいてあることはあまり変わらないな。提案手法と比較したのは最近の RNN ベースの手法と、matrix factorization 手法と、伝統的な手法だとはかいてある。あと提案手法の利点は季節構造を明示的にモデリングできるとこらしい。…あのさ、この論文ちょいちょい Appendix を参照ってあるんだけど、Appendix なくない?

f:id:cookie-box:20180305232608p:plain:w60

確かにないですね…あ、Proceedings のページの "PDF" ではなくて "Supplemental" というリンク先のファイルには Appendix が付いていますよ。

f:id:cookie-box:20180305231302p:plain:w60

マジで? よかった、メールで問い合わせないといけないかと思った。だったら俺英作文できないからハルナに頼まなきゃだし。

f:id:cookie-box:20180305232608p:plain:w60

春名さんは英語を話せるけど書けないんじゃ…というか春名さんにまず日本語で趣旨が伝えられるでしょうか。「論文の Appendix が見当たらないのですが」という質問をしたいんですよね?

f:id:cookie-box:20180305231302p:plain:w60

駄目だその日本語ハルナに通じないな…。まあ Appendix あったからいいや。Intoroduction を読むと時系列予測でも特に需要予測は現実で重要だってかいてるな。確かに "electricity" と "traffic" も需要予測っぽいな。道路の需要っていうとちょっと違和感あるけど。とにかくそういう予測では正確で最新のデータに基づいた予測値がほしいと。それで、時系列のトレンドや季節性を捉えるモデルの最たる例が状態空間モデルだけど、状態空間モデルは時系列の構造があらかじめ結構わかっているときには有利で、解釈もしやすいし、データ効率?のいい学習もできるけど、一方でじゅうぶんな長さの学習データが必要だし、あと巨大で複雑な時系列データを扱うには説明変数を決めるのが大変すぎだって。さらに、伝統的な状態空間モデルでは時系列をばらばらにフィッティングするから類似の時系列が共有するパターン?を推論できない? そうなの?

f:id:cookie-box:20180305232608p:plain:w60

状態空間モデル自体はもちろん多変量時系列を扱うことができますし、その多変量の分布をトラッキングできますが…「類似の時系列」というのがポイントなんでしょうか。「類似の時系列」というのは明らかに「いろいろな顧客の消費電力の時間変化のパターン」「いろいろな道路の混み具合の時間変化のパターン」のようなものを指していますよね。このような時系列を予測するとき、「人が活動しない夜中は消費電力が少ないだろう」などという特徴を考慮してシステムモデルをつくるんだと思います。でもそのシステムモデルを370人の顧客それぞれに適用したのでは、370人分のデータがある意味があまりないんだと思います。もちろん「顧客Aと顧客Bにの消費電力には強い正の相関がある」などというのは結果的にトラッキングされますけど、「顧客Aと顧客Bにの消費電力が強い正の相関をもつのはどのようなときである」というのを積極的に学習するわけではない、というようなイメージでしょうか…。

f:id:cookie-box:20180305231302p:plain:w60

ふーん? それで逆に深層ニューラルネットワークは高次元の特徴を抽出することができて、あまり人が手間をかけることなく時系列を横断する複雑なパターンを特定できるって。本当かなー。ただこの手のモデルはあまり時系列の構造に何か想定をおいてるものじゃないから、正確な学習データを得るにはやっぱり長期間分の学習データが必要で、解釈もしにくいし、平滑化みたいなことをしたくても織り込みにくいって。

f:id:cookie-box:20180305232608p:plain:w60

「伝統的なモデルは解釈性と可塑性に優れるが手間を要し表現力に欠け、ニューラルネットは手間がかからず表現力に長けるが解釈性と可塑性に欠ける」という導入に幾度とない既視感がある気がしますね。

f:id:cookie-box:20180305231302p:plain:w60

それで便利なモデルを考えたなら別にいいだろ! でこの論文はもちろん状態空間モデルとニューラルネット、とりわけ RNN みたいだけど、を組み合わせるんだけど、「線形状態空間モデルを RNN でパラメタライズする」らしい…これだけじゃ意味わかんないな。でも、個別の時系列は状態空間モデルでモデリングして、時系列間を横断する構造を RNN でモデリングするっぽい?

f:id:cookie-box:20180305232608p:plain:w60

人工データで検証しているということですから、どのような人工データを生成しているかで著者らがどのような構造の学習を意図しているのかわかるでしょうね。

f:id:cookie-box:20180305231302p:plain:w60

提案モデルは状態空間モデルをつかってるから解釈もしやすいし、状態空間モデルに入れるモデルの想定によって少ないデータにも大量のデータにも適応できる。ってデータが少なくても多くてもいいって状態空間モデルパート任せだったのか…。まあだから過学習を防げるって。データが少ないと過学習しやすいから、データが少ないなりの想定を状態空間モデルに入れとけってことだよな…。とまあイントロダクションとしてはこんな感じで、2節は先行研究だけど。

f:id:cookie-box:20180305232608p:plain:w60

短いので確認しておきましょうか。Hyndman et al. や Durbin and Koopman が状態空間モデルについてまとめており、最近では機械学習分野で状態空間モデルに関する研究もみられると。状態空間モデルと RNN の組合せというのはこれまでにも提案されていて、以下の片方もしくは両方に該当するということです。

  • ガウシアンではない複雑な尤度関数に拡張する。
  • MLP によって時間発展に非線形性をもたせる、または状態空間モデルと RNN によって特定した時間変化する遷移行列を組み合わせる。
というか状態空間モデルと深層学習の組合せの例がかなり挙げられていますが、論文上の記述からのみではアイデアがわかりづらいのでタイトルだけ列挙しておきます。
  • Deep Markov Model (DMM) [Krishnan2017, Krishnan2015]
  • Stochastic RNNs [Fraccaro2016]
  • Variational RNNs [Chung2015]
  • Latent LSTM Allocation [Zaheer2017]
  • State-Space LSTM [Zheng2017]
  • 状態空間モデルの教師なし学習 [Karl2017]
  • Kalman Variational Auto-Encoder (KVAE) [Fraccaro2017] ― 本論文での提案手法に最も近い。
  • 上と似たアイデア [Johnson2016]
この下から2番目のモデルが本論文の提案手法と似ているらしいですが、こちらの手法では RNN で固定された K 個のパラメータを学習しているのに対して、本論文の手法では RNN で直接的に状態空間モデルのパラメータを学習するのでハイパーパラメータチューニングが不要とありますが…よくわからないので先に進みましょう。

f:id:cookie-box:20180305231302p:plain:w60

3節は一般的な状態空間モデルによる予測の話らしいけど、表記や前提をつかむために読んでおいた方がいいな。いま  \{ z_{1:T_i}^{(i)}\}_{i=1}^N っていう N 本の時系列があるっぽい。 \{ z_{1:T_1}^{(1)} \}, \; \{ z_{1:T_2}^{(2)} \}, \cdots, \; \{ z_{1:T_N}^{(N)}\} ってことだな。それで  z_t^{(i)} \in \mathbb{R} だな。脚注に、各時系列は等間隔だけど1本目の時系列の時刻 t=1 と2本目の時系列の時刻 t=1 はそろってなくてもいいってあるけど、まあ気にしなくていいかな? それで、それに伴う説明変数もあって  \{ x_{1:T_i + \tau}^{(i)}\}_{i=1}^N こう。こっちの  x_t^{(i)} はベクトルで  x_t^{(i)} \in \mathbb{R}^D だ。なんで  x_t^{(i)} の方が  z_t^{(i)} より \tau だけ長く手元にあるんだって思ったけど、\tau ステップ先が予測対象で、 x_t^{(i)} はその \tau ステップ先までの値が手に入っていて、その状態で  z_{T_i + 1:T_i + \tau}^{(i)} を予測するのか。以下の確率分布を求めたいらしい。

p \left( z_{T_i + 1:T_i + \tau}^{(i)}  \, \middle| \, z_{1:T_i}^{(i)} , \, x_{1:T_i + \tau}^{(i)}   ; \, \Phi \right)
縦棒の右側に他の時系列は入ってこないんだな。でもモデルのパラメータ \Phi は他の時系列とシェアされてるらしい。それで仮定として、  x_{1:T_i}^{(i)}\Phi が given の下では N 本の時系列  \{ z_{1:T_i}^{(i)}\}_{i=1}^N は時系列間で独立とするって。この仮定を置いちゃうと時系列間の相関構造はモデリングできないけど…って、あ!!

f:id:cookie-box:20180305232608p:plain:w60

どうしたんですか?

f:id:cookie-box:20180305231302p:plain:w60

時系列間が独立と仮定しちゃうと相関構造をモデリングはできないけど、パラメータ \Phi が共通だから、時系列間で「統計的な強度」をシェアできないわけじゃない、ってかいてあるんだよ。statistical strength って、PRML の6章の最後で出てきたじゃん!

このときも、ニューラルネットの各出力変数は隠れ層を共有するから互いに「統計的な強度」を借りることができるっていってた。

f:id:cookie-box:20180305232608p:plain:w60

多変量を予測するモデルは、個別にモデリングするよりもまとめてモデリングして共通因子をもたせる方が個々の変数をより確かに予測できるという意味なんでしょうね。

f:id:cookie-box:20180305231302p:plain:w60

それで、状態空間モデルは典型的には1本1本の時系列に対して適用するって。こういう電力消費量の予測みたいな文脈ではってこと? だって気象予報の状態空間モデルでは全然そんなことないよな…。まあともかく、こういうシステムモデルを考える。

 l_t = F_t l_{t-1} + g_t \varepsilon_t , \qquad \varepsilon_t \sim N(0,1)
l_t が状態で F_t が遷移行列で、ノイズもガウシアンだし、カルマンフィルタが適用できる状況だな。それで観測モデルがこう。やっぱり線形でガウシアンな観測を考えてる。
 z_t = y_t + \sigma_t \epsilon_t, \quad y_t = a_t^{\top} l_{t-1} + b_t , \quad \epsilon_t \sim N(0,1)
それで、特定しなきゃいけないパラメータは  \Theta_t = (\mu_0, \Sigma_0, F_t, g_t, a_t, b_t, \sigma_t) で、これを特定するには例えば周辺尤度を最大化するようなパラメータを選ぶ。周辺尤度は以下だな。
 \displaystyle p_{SS}(z_{1:T}| \Theta_{1:T}) \equiv p(z_1|\Theta_1) \prod_{t=2}^{T} p(z_t | z_{1:t-1}, \Theta_{1:t})
それでパラメータが特定できたら z_{1:T} の分布は l_{1:T}積分すれば出るけど、今回の条件では z_{1:T} の分布は l_{1:T}積分しなくてもカルマンフィルタで解析的に求まるな。でも、せっかく N 本の時系列があるのに1本1本予測してたんじゃ時系列間で情報が共有されないから、学習データが少なかったりノイズが大きかったりしても適用できるようにこれから RNN と組み合わせるって。

f:id:cookie-box:20180305232608p:plain:w60

それで4節がその Deep State Space Models…各時刻に  x_t^{(i)} を入力し、各時刻の出力が  \Theta_t^{(i)} となるように RNN を学習するということですか。論文タイトルが Deep State Space Model"s" と複数形なのも N 本の各時系列の状態空間モデルを同時に学習するからだったんですね。ネットワークの学習の損失関数は周辺尤度ですね。RNN は時系列からその時刻までの特徴を抽出するのに適したモデルなので、学習対象を状態空間モデルそのものにするというのは理解できる気がしますが…これ、"electricity" と "traffic" における「説明変数(状態)」とは何なんでしょう?

f:id:cookie-box:20180305231302p:plain:w60

4節端折りすぎだろ! この論文のメインパートなのに! 5節の検証パートは Appendix を合わせてみないとわからないなこれ…。これたぶん"electricity" と "traffic" ではどっちも状態空間モデルでも季節性に基づくモデルを採用してて、トラッキングしてる状態は 31 次元で、そのうち 24 次元は1日の中での時間の季節成分、7 次元は曜日の季節成分をあらわしてる。詳細はこっちの論文みろって。この DeepAR は比較対象にしてる手法の1つだな。

次回があればつづく

electricity のプロット
import numpy as np
import pandas as pd
import datetime

df = pd.read_csv('LD2011_2014.txt/LD2011_2014.txt', sep=';', dtype='object', nrows=600)
df.columns.values[0] = 'time' # 時刻カラムだけ列名がないので付ける

# 時刻カラムの処理
df['time'] = df['time'].map(lambda x: datetime.datetime.strptime(x, '%Y-%m-%d %H:%M:%S'))

# 時刻でないカラムの処理 (このデータは小数点がカンマになっている)
for i in range(1, 371):
    df.iloc[:,i] = df.iloc[:,i].map(lambda x: x.replace(',', '.')).astype(np.float64)
    #if df.iloc[:,i].sum() > 0.0:
    #    print(i)

%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 14, 6
rcParams['font.size'] = 16

plt.plot(df['time'], df['MT_191'])
plt.plot(df['time'], df['MT_192'])
plt.plot(df['time'], df['MT_193'])
plt.plot(df['time'], df['MT_194'])
plt.plot(df['time'], df['MT_195'])
plt.legend()
plt.xlabel('time')
plt.ylabel('electricity')
plt.show()
traffic のプロット
import numpy as np

x = []
f = open('PEMS-SF/PEMS_train')
line = f.readline()
n = 0
while line:
    rows = line.strip('[]\n').split(';') # ファイルの各行は大カッコでくくられ、行の境目がセミコロンで、行内が空白で区切られた行列
    x_ = []
    for row in rows:
        x_.append(np.array(row.split(' ')).astype(np.float64))
    x.append(x_)
    n = n + 1
    if n == 11:
        break
    line = f.readline()
f.close()

x = np.array(x)
#print(x)
#print(x.shape)

%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 14, 6
rcParams['font.size'] = 16

# 0番目が水曜日で10番目が日曜日( PEMS_trainlabels をみる)
plt.plot(np.arange(0, 24, step=1/6), x[0][0], label='Wednesday Sensor0', color='dodgerblue')
plt.plot(np.arange(0, 24, step=1/6), x[0][1], label='Wednesday Sensor1', color='dodgerblue', linestyle='dashed')
plt.plot(np.arange(0, 24, step=1/6), x[0][2], label='Wednesday Sensor2', color='dodgerblue', linestyle='dotted')
plt.plot(np.arange(0, 24, step=1/6), x[10][0], label='Sunday Sensor0', color='darkorange')
plt.plot(np.arange(0, 24, step=1/6), x[10][1], label='Sunday Sensor1', color='darkorange', linestyle='dashed')
plt.plot(np.arange(0, 24, step=1/6), x[10][2], label='Sunday Sensor2', color='darkorange', linestyle='dotted')
plt.legend()
plt.xlabel('hour')
plt.ylabel('occupancy rate')
plt.show()