graphvizで決定木を可視化

前回の記事でgraphvizをインストールしたので、早速決定木を可視化してみましょう。
サンプルなので、モデル自体は適当に作ります。(データもいつものirisです)


# ライブラリインポート
import graphviz
from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier
from sklearn.tree import export_graphviz
# データの読み込みとモデルの学習
iris = load_iris()
X = iris.data
y = iris.target
clf = DecisionTreeClassifier(min_samples_split=5)
clf.fit(X, y)

これで、モデルができたので可視化してみましょう。
参考にするコードはscikit-learnのサンプルがいいと思います。
ここではblogに挿入することを考慮して、pdfではなくpng画像に出力しました。


dot_data = export_graphviz(
                        clf,
                        class_names=iris.target_names,
                        feature_names=iris.feature_names,
                        filled=True,
                        rounded=True,
                        out_file=None
                    )
graph = graphviz.Source(dot_data)
graph.render("iris-tree", format="png")

結果がこちらです。
わかりやすく可視化できましたね。

macにgraphvizをインストールする

決定木の可視化などに使うgraphvizをMacにインストールする手順のメモです。
ついでにpythonライブラリもインストールします。

まず、本体の方はドキュメントにある通り、homebrewでインストールできます。


brew install graphviz

pythonライブラリの方はこちら。 pipコマンドでインストールします。


pip install graphviz

動作確認をかねて Quickstartのコードを動かしてみましょう。


>>> from graphviz import Digraph
>>> dot = Digraph(comment='The Round Table')
>>> dot  #doctest: +ELLIPSIS
>>> dot.node('A', 'King Arthur')
>>> dot.node('B', 'Sir Bedevere the Wise')
>>> dot.node('L', 'Sir Lancelot the Brave')
>>> dot.edges(['AB', 'AL'])
>>> dot.edge('B', 'L', constraint='false')
>>> print(dot.source)  # doctest: +NORMALIZE_WHITESPACE
>>> dot.render('test-output/round-table.gv', view=True)

これでpdfファイルを開いた時に3個の丸を含むグラフが表示されたらOKです。

(2020/3/21 追記)
このページへのアクセスが非常に多い理由が不明だったのですが、
conda でインストールして、うまく動かなかった人がいるのではないかと思ったので補足します。

condaの方にも、 graphviz と言うパッケージがあり、
conda install graphviz でインストールできますが、
これはgraphvizの本体が含まれないラッパーのみで、
しかも本体を別途 brew でインストールしてもパスが通らないことがあるようです。

conda で導入する場合は
conda install python-graphviz
を使います。

ほかのサイトでは、
brew で本体を入れた後に, graphviz と python-graphviz をそれぞれ入れている例を見かけますが、
試したところ、brewも無しで、
conda install python-graphviz
だけ、行えば動きました。

kerasのモデルの中間層の出力を可視化してみる

ディープラーニングのモデルを作成したとき、中間層の出力が気になることがよくあります。
きちんと活性化しているかとか、相関が高すぎて意味がないユニットが多くないかとか、
どんな条件の時に活性するのかなど、確認したい内容は時により様々です。

kerasの場合、学習済みのモデルの層を取り出して新しいモデルを作成することで中間層の出力を確認できます。

中間レイヤーの出力を得るには?

試しに以前下記の記事で作ったモデルでやってみましょう。
CNNで手書き数字文字の分類

公式ドキュメントに紹介されていたのと少し違う方法ですが、普通にSequentialモデルに学習済みの層を一個追加したら動いたので、
その方法で行います。
一層目には16ユニットあるのですが、そのうち2このユニットについて、出力を可視化しました。

# 学習済みモデルの1層目だけ取得してモデルを作成する


model_2 = Sequential()
model_2.add(model.layers[0])

# 元画像と1層目の出力2個を可視化
fig = plt.figure(figsize=(18, 30))
for i in range(5):
    # print(y_test[i].argmax())
    ax = fig.add_subplot(6, 3, 3*i+1)
    ax.imshow(X_test[i][:, :, 0], cmap='gray_r')
    ax = fig.add_subplot(6, 3, 3*i+2)
    ax.imshow(model_2.predict(X_test[i:i+1])[0][:, :, 0], cmap='gray_r')
    ax = fig.add_subplot(6, 3, 3*i+3)
    ax.imshow(model_2.predict(X_test[i:i+1])[0][:, :, 1], cmap='gray_r')
plt.show()

出力がこちらです。

真ん中の列の出力は横線の下辺に反応していることや、右側の列の結果は中抜き文字のような形で反応しているのがわかりますね。

ちなみに、それぞれのユニットのウェイト(バイアスは除く)を可視化すると次のようになります


fig = plt.figure(figsize=(5,10))
for i in range(2):
    w = model_2.get_weights()[0][:, :, 0, i].reshape(3, 3)
    ax = fig.add_subplot(2, 1, i+1)
    ax.imshow(w, cmap='gray_r')
plt.show()

イメージした通りのウェイトでした。

PrestoのURL関数

Webサービスのアクセスログデータを集計するとき、URLを扱う場面はよくあります。
WHERE句で特定の条件のURLを絞り込んだり、
SELECT句でパラメーターやパスごとに数えたりする場面はよくありますが、
Prestoにはそのようなときに便利な関数が準備されていいます。

ドキュメント 6.18. URL Functions

具体的には下記の関数群です。

  • url_extract_fragment(url)
  • url_extract_host(url)
  • url_extract_parameter(url, name)
  • url_extract_path(url)
  • url_extract_port(url)
  • url_extract_protocol(url)
  • url_extract_query(url)
  • URLを次の形式とすると、それぞれの関数で url_extract_xxxx の xxxx 部分を文字列として返してくれます。
    [protocol:][//host[:port]][path][?query][#fragment]

    これを知るまでは LIKE や正規表現を使うことが多く、稀に予期せぬパターンがマッチしてしまったりしていて困っていたのですが、
    これらのURL関数群を使うようになってからは、そのようなエラーは起こりにくくなりましたし、
    クエリも意図がわかりやすいものになりとても助かっています。

PrestoのCOALESCE関数でNULL値を別の値に置き換える

Prestoの便利な関数の紹介です。
Prestoに限らず、SQLデータを抽出するときにNULL値を別の値に置き換えたいケースは頻繁に発生します。
例えば、NULLならば0、その他の値はそのまま出力したいときは下記のような書き方になります。


CASE
    WHEN column_name1 IS NULL THEN 0
    ELSE column_name1
END AS column_name1

これを、下記のように少しだけスッキリ書くことがでます。


COALESCE(column_name1, 0) AS column_name1

ドキュメントはこちら

この他、次のように引数を複数並べて書くと、左から順番に評価して最初にNULLでなかった列の値を取得できます。
完全外部結合するときなどに便利ですね。


COALESCE(column_name1, column_name2, column_name3, 0) AS column_name1

pythonで編集距離(レーベンシュタイン距離)を求める

ごく稀にですが、文字列同士の編集距離を求める必要が発生するのでその時のメモです。

編集距離(レーベンシュタイン距離)とは、二つの文字列がどの程度異なっているかを表す距離の一種です。
Wikipediaにも解説があります。

一方の文字列に対して、1文字の挿入、削除、置換を最低何回施せばもう一方の文字列に等しくなるかで定まります。

pythonでこれを求めるときは、python-Levenshtein というライブラリが使えます。

インストール


pip install python-Levenshtein

使い方


>>> import Levenshtein
>>> text1 = 'Levenshtein'
>>> text2 = 'Lenvinsten'
>>> Levenshtein.distance(text1, text2)
4

CNNで手書き数字文字の分類

以前の記事で読み込んだ手書き数字文字データ(MINIST)を使って、0~9の数字を判定するモデルを作ってみます。

kerasのサンプルコードもあるのですが、せっかくなので少しだけパラメーターなどを変えてやってみましょう。

最初にライブラリをインポートしてデータを準備します。
前処理として配列の形をConv2Dのinputに合わせるのと、
0〜1への正規化、 ラベルの1hot化を行います。


# ライブラリの読み込み
from keras.datasets import mnist
from keras.utils import to_categorical
from keras.models import Sequential
from keras.layers import Dense, Dropout, Flatten
from keras.layers import Conv2D, MaxPooling2D
from keras.callbacks import EarlyStopping
import matplotlib.pyplot as plt
from sklearn.metrics import classification_report

# データの読み込み
(data_train, target_train), (data_test, target_test) = mnist.load_data()

# Conv2D の inputに合わせて変形
X_train = data_train.reshape(-1, 28, 28, 1)
X_test = data_test.reshape(-1, 28, 28, 1)

# 特徴量を0~1に正規化する
X_train = X_train / 255
X_test = X_test / 255

# ラベルを1 hot 表現に変換
y_train = to_categorical(target_train, 10)
y_test = to_categorical(target_test, 10)

続いてモデルの構築です


# モデルの構築
model = Sequential()
model.add(Conv2D(16, kernel_size=(3, 3),
                 activation='relu',
                 input_shape=(28, 28, 1)))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))
model.add(Conv2D(32, (3, 3), activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Dropout(0.25))
model.add(Flatten())
model.add(Dense(64, activation='relu'))
model.add(Dropout(0.5))
model.add(Dense(10, activation='softmax'))
model.compile(
    loss="categorical_crossentropy",
    optimizer="adam",
    metrics=['accuracy']
)
print(model.summary())

# 以下、出力
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
conv2d_1 (Conv2D)            (None, 26, 26, 16)        160       
_________________________________________________________________
max_pooling2d_1 (MaxPooling2 (None, 13, 13, 16)        0         
_________________________________________________________________
dropout_1 (Dropout)          (None, 13, 13, 16)        0         
_________________________________________________________________
conv2d_2 (Conv2D)            (None, 11, 11, 32)        4640      
_________________________________________________________________
max_pooling2d_2 (MaxPooling2 (None, 5, 5, 32)          0         
_________________________________________________________________
dropout_2 (Dropout)          (None, 5, 5, 32)          0         
_________________________________________________________________
flatten_1 (Flatten)          (None, 800)               0         
_________________________________________________________________
dense_1 (Dense)              (None, 64)                51264     
_________________________________________________________________
dropout_3 (Dropout)          (None, 64)                0         
_________________________________________________________________
dense_2 (Dense)              (None, 10)                650       
=================================================================
Total params: 56,714
Trainable params: 56,714
Non-trainable params: 0
_________________________________________________________________

そして学習です。


# 学習
early_stopping = EarlyStopping(
                        monitor='val_loss',
                        min_delta=0.0,
                        # patience=2,
                )

history = model.fit(X_train, y_train,
                    batch_size=128,
                    epochs=30,
                    verbose=2,
                    validation_data=(X_test, y_test),
                    callbacks=[early_stopping],
                    )

# 以下出力
Train on 60000 samples, validate on 10000 samples
Epoch 1/30
 - 18s - loss: 0.6387 - acc: 0.7918 - val_loss: 0.1158 - val_acc: 0.9651
Epoch 2/30
 - 18s - loss: 0.2342 - acc: 0.9294 - val_loss: 0.0727 - val_acc: 0.9772
Epoch 3/30
 - 17s - loss: 0.1827 - acc: 0.9464 - val_loss: 0.0571 - val_acc: 0.9815
Epoch 4/30
 - 18s - loss: 0.1541 - acc: 0.9552 - val_loss: 0.0519 - val_acc: 0.9826
Epoch 5/30
 - 18s - loss: 0.1359 - acc: 0.9598 - val_loss: 0.0420 - val_acc: 0.9862
Epoch 6/30
 - 17s - loss: 0.1260 - acc: 0.9620 - val_loss: 0.0392 - val_acc: 0.9880
Epoch 7/30
 - 18s - loss: 0.1157 - acc: 0.9657 - val_loss: 0.0381 - val_acc: 0.9885
Epoch 8/30
 - 19s - loss: 0.1106 - acc: 0.9673 - val_loss: 0.0349 - val_acc: 0.9889
Epoch 9/30
 - 17s - loss: 0.1035 - acc: 0.9694 - val_loss: 0.0359 - val_acc: 0.9885

学習の進み方をプロットしておきましょう。


# Epoch ごとの正解率と損失関数のプロット
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(2, 1, 1, title="loss")
ax.plot(history.epoch, history.history["loss"], label="train_loss")
ax.plot(history.epoch, history.history["val_loss"], linestyle="-.", label="val_loss")
ax.legend()
ax = fig.add_subplot(2, 1, 2, title="acc")
ax.plot(history.epoch, history.history["acc"], label="train_acc")
ax.plot(history.epoch, history.history["val_acc"], linestyle="-.", label="val_acc")
ax.legend()
plt.show()

val_acc は結構高い値を出していますが、一応クラスごとの成績も評価しておきましょう。


# 評価
y_predict = model.predict_classes(X_train)
print(classification_report(target_train, y_predict))

# 以下出力
             precision    recall  f1-score   support

          0       0.99      1.00      0.99      5923
          1       0.99      1.00      0.99      6742
          2       0.99      0.99      0.99      5958
          3       0.99      0.99      0.99      6131
          4       0.99      0.99      0.99      5842
          5       0.99      0.99      0.99      5421
          6       0.99      0.99      0.99      5918
          7       0.98      0.99      0.99      6265
          8       0.99      0.97      0.98      5851
          9       0.98      0.99      0.98      5949

avg / total       0.99      0.99      0.99     60000

ほとんど適当に作ったモデルでしたが、
ほぼほぼ正解できていますね。

pickleを使ってpythonのオブジェクトをファイルに保存する

(注)この記事はscikit-learnのモデルをファイルに保存することを念頭に書いていますが、
pickle自体はscikit-learnのモデル以外のものも直列化してファイルに書き出すことができるモジュールです。

以前の記事で、kerasで作成したmodelを保存したり読み込んだりする方法を書きました。
今回はscikit-learnで作ったモデルを保存してみます。
kerasには専用の関数が用意されていたのですが、scikit-learnにはありません。
そのため、他の方法が必要です。
そこでpython標準ライブラリの pickleが使えます。
ドキュメント

利用方法は、ドキュメントのpickle.dumppickle.loadの説明と、一番下の使用例が参考になります。

clfという変数に、学習済みのモデルが格納されているという想定で、保存と読み込みのコード例を紹介します。
また、保存するファイル名は何でも良いのですが、サンプルコードではclf.pickleとします。

まずは保存。


import pickle
with open("clf.pickle", "wb") as f:
    pickle.dump(clf, f)

次に読み込み。


import pickle
with open("clf.pickle", "rb") as f:
    clf = pickle.load(f)

これで、一度学習したモデルを読み込んで、予測に活用することができます。
scikit-learnで学習したモデルを本番運用するならばほぼ必須の技術です。
(pickle以外の方法を使うという手もありますが、何らかの形での保存と読み込みの手段が必要です)

ブログ記事数が50を超えていたので振り返り

このブログをはじめて1ヶ月半ほどたち、記事数はいつの間にか50記事を超えていました。
想定外だった気づきや想定外だったこともあるので、この辺で一度振り返りをやっておこうと思います。
とりあえずこのブログについて感じていることをそのまま書き出しておきます。

想像以上に勉強になる

元々は、自宅PC、職場PC、個人アカウントのサーバー、職場のサーバーなどに散らばっていたメモ書きを
一箇所に集約して参照できたら便利だということと、
それがネット上にあれば自分より後に勉強を始めた人の役に立つだろうという二つの目的でこのブログを始めました。
要は、最初は自分が知ってることとか調べ終わったことを書き写すだけで新しい学びが増えることはあまり期待していましせんでした。
しかし、改めて記事を書くとなるとドキュメントを読み直したり、最新バージョンの環境でテスト実行したりするので、
新たな気づきが多くあり自分の勉強になりました。
numpyやpandasなんて2年弱使っているのですが、
これまでは、本や誰かが書いたサンプルコードばかり読んで、公式ドキュメントをまともに見ていませんでした。
やはり公式のドキュメントを読むとより効率のいい使い方がわかったり、新しく知る機能があったりと非常に勉強になります。

記事を書くのは思ったより大変

上の内容とかぶりますが元々は既存のメモ書きを元ネタに記事を書こうとしていたので、もっとサクサクと記事を増やせると思っていました。
しかし、記事用に調べ直したり、できるだけ綺麗なスタイルでコードを書いたり、画像を準備したりすると、
既知の内容でも結構時間がかかります。新しく知ったテーマならなおさらです。
今までいろんな人のブログや記事、投稿を読んできましたが、それらを書いてくださっていた
皆さんへの感謝の気持ちと尊敬の念が高まりました。

想定していた内容の記事がまだ書けてない

データサイエンティストになってからユーザーの行動分析や自然言語処理に多く取り組んできたので、
その分野でもっと高度な記事を量産できると思っていました。
実際、書こうと思っているネタは色々ありますので、そのうち書きます。

しかし実際は最初の方の記事はAWSにwordpressを構築するために色々調べた内容ばかりになり、
その後もそれこそデータサイエンティストになったばかりの頃(というよりもpythonを触り始めた頃)のメモ書きを見返して、
非常に基本的な内容の記事ばかり書いています。
ただ、前述の通りそれが勉強になっているので当分はそのまま続けようと思っています。
基本的な内容でもきっと誰かの役に立つこともあるでしょうし。

まだアクセスが全然増えない

最近は、1日2〜3人くらいは自分以外の人が来てくれているようなのですが、
そのレベルで留まって全然アクセスが伸びていません。
昔やっていたyahooやfc2のブログではもっと早い段階で訪問者が増えたのですが、
それぞれのサービスにあった新着記事のフィードに乗る効果やSEOのおかげだったのでしょう。

今の段階で訪問が増えても、基本的な内容しか書いてなくて若干恥ずかしいというのもあるので、
あんまり積極的なアクセスアップ施策は取らず、当分地道に記事を増やしていきたいと思っています。

今後の計画

とりあえず今のところは1日1記事のペースを維持できていますが、
より高度な内容を書こうと思うとこのペースを維持するのは厳しいとも感じています。
まだまだネタはあるので、当分1日1記事ペースを頑張りたいですがそれをずっと続けるのはおそらく無理です。
どこかのタイミングで、ちゃんと時間かけて独自に調べた内容を週1回くらいのペースで
上げていくスタイルに移行できた方が読む人のためになるのかなとも思っています。

pandasでgroupbyした時に複数の集計関数を同時に適用する

前の記事の続きです。
pandasでデータフレームをgroupbyした時に使える集計関数
ドキュメントのこの記事で参照した部分のすぐ下に、
Applying multiple functions at once
という段落があります。
実はこれ初めて知りました。
今までグルプごとに個数と、平均と、標準偏差を計算したい、みたいな時は、
groupbyして集計を個別に実施して、その結果をmergeするという非常に面倒なことをずっとやっていました。

それが、aggというのを使うと一発でできるようです。


import pandas as pd
from sklearn.datasets import load_iris

# データフレームの準備
iris = load_iris()
df = pd.DataFrame(iris.data, columns=iris.feature_names)
df["target"] = iris.target
df["target_name"] = df.target.apply(lambda x:iris.target_names[x])
del df["target"]

df.groupby("target_name").agg(["count", "mean", "std"])

出力されるのが次です。(ブログのレイアウトの都合上画像で貼り付けます。)

これは便利です。
また、DataFrameのカラム名が2段になっています。
これをみて、indexだけではなく実はcolumnsでも、MultiIndexが使えることを知りました。