カテゴリー: 機械学習

ライブラリやアルゴリズムなど、機械学習に関すること

投稿日:

scikit-learnでK-分割交差検証

普段、機械学習のパラメーターを最適化するために
K分割交差検証をするときはGridSearchCVで済ましているのですが、
今回別の目的があって、K-分割交差検証のような分割が必要になりました。
参考:scikit-learn でグリッドサーチ

そこで、GridSearchCVの中でも使われているKFoldを使ってみたのでその記録です。
ドキュメントはこちら。
sklearn.model_selection.KFold

注意点は、splitした時に得られる結果が、データそのものではなくデータのインデックスである点くらいです。

実行してみるために、まずダミーデータを準備します。


from sklearn.model_selection import KFold
import numpy as np

# ダミーデータの準備
X = np.random.randint(20, size=(10, 2))
y = np.random.choice([0, 1], size=10)
print(X)
'''
[[ 2 12]
 [ 6  3]
 [ 0  3]
 [15  9]
 [19  0]
 [12 14]
 [ 0  0]
 [12 18]
 [ 1  9]
 [ 3 14]]
'''
print(y)
# [0 0 0 1 1 1 0 1 0 0]

次に、 KFoldを動かしてみましょう。
まず挙動を確認したいので、splitして得られるイテレータの中身を表示してみます。


kf = KFold(5, shuffle=True)
for train_index, test_index in kf.split(X):
    print("train_index:", train_index)
    print("test_index:", test_index)
    print()  # 1行開ける

# 以下出力
'''
train_index: [1 2 3 4 5 6 8 9]
test_index: [0 7]

train_index: [0 1 2 3 4 6 7 9]
test_index: [5 8]

train_index: [0 1 3 4 5 6 7 8]
test_index: [2 9]

train_index: [0 1 2 5 6 7 8 9]
test_index: [3 4]

train_index: [0 2 3 4 5 7 8 9]
test_index: [1 6]

'''

ご覧の通り、0〜9の値を、5つのグループに分けて、順番に一個をテスト用にしながら交差検証用のデータセットのインデックスを返してくれています。
インデックスを元にデータを分割するには、つぎのようにして別の変数に取り出すと以降のコードが読みやすくなります。


X_train, X_test = X[train_index], X[test_index]
y_train, y_test = y[train_index], y[test_index]
投稿日:

LightGBMを動かしてみる

今回はほぼ動作確認だけです。
(前回の記事の最後に書いた通り、僕の環境ではimport 時にWarningが出るのでそもそも動くかどうかが不安だったので。)

LightGBMには scikit-learn形式のAPIが実装されていて、ほぼ同じように使う事もできます。
Scikit-learn API

とりあえずbostonのデータセットを使って、回帰問題を解いてみましょう。
動作確認なのでパラメーターはデフォルトです。
最後の評価は平方二乗誤差で行なってます。


import lightgbm
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

# データの準備
boston = load_boston()
X = boston.data
y = boston.target
X_train, X_test, y_train, y_test = train_test_split(
                                    X, y, test_size=0.2, random_state=0
                                )
# デフォルトパラメーターでオブジェクト生成
lgbm = lightgbm.LGBMRegressor()
# 学習
lgbm.fit(X_train, y_train)
# 評価
y_predict = lgbm.predict(X_test)
print(mean_squared_error(y_test, y_predict))
# 24.34062206680735

無事に動いたようです。

投稿日:

pythonでトピックモデル(LDA)

この前の記事で、scikit-learnのニュース記事のジャンルをロジスティック回帰で予測するというモデルを作ってみました。

参考:scikit-learnのニュースデータをロジスティック回帰で分類

今回はアプローチを変えて、トピックモデルを試してみようと思います。
どちらかというと、20newsのデータセットでもう少し何かやりたいというのが主目的で、
トピックモデルの理論的な説明については今回は省略します。
興味のあるかたへは、講談社から出ている岩田具治先生の、 トピックモデル (機械学習プロフェッショナルシリーズ)
が非常にわかりやすかったのでおすすめです。ページ数も少なめでありがたい。
(数式が多くて書くのが大変なのですがゆくゆくは時系列分析みたいにこのブログでも説明したい。)

さて、pythonでトピックモデルを実装するには gensim を使うのが一般的のようです。
gensim topic modelling for humans
ただ、今回はいつも使っているscikit-learnでやってみました。
(gensimはword2vec等で使ってるのですがscikit-learnに比べると少し苦手。)

scikit-learnでトピックモデルを実装するために読むドキュメントはこちら。
sklearn.decomposition.LatentDirichletAllocation
Topic extraction with Non-negative Matrix Factorization and Latent Dirichlet Allocation

サンプルコードと同じことをしてもしょうがないので、少し工夫をしています。
・サンプルデータのカテゴリーを前回の記事同様に5個に絞る(その代わりそのカテゴリの全データを使用)
・カテゴリーごとに各文章のトピッックを可視化

前置きが長くなりましたが、やってみましょう。
必要ライブラリーのインポートとデータの読み込み


from sklearn.feature_extraction.text import CountVectorizer
from sklearn.decomposition import LatentDirichletAllocation
from sklearn.datasets import fetch_20newsgroups
import matplotlib.pyplot as plt
import numpy as np

remove = ('headers', 'footers', 'quotes')
categorys = [
        "rec.sport.hockey",
        "soc.religion.christian",
        "sci.med",
        "comp.windows.x",
        "talk.politics.mideast",
    ]
twenty_news = fetch_20newsgroups(
                                subset='all',
                                remove=remove,
                                categories=categorys
                            )
X = twenty_news.data

続いて、単語の出現頻度を数え、LDAのモデルを構築して学習します。
トピック数は カテゴリー数と同じ5でも試したのですが、
どうやら6か7にして、あまり重要でない単語を引き受けるトピックを作った方が納得性の高いものになりました。
サンプルコードは20カテゴリーを10トピックでうまく処理できているのに何故だろう？


# 単語の出現頻度データを作成
tf_vectorizer = CountVectorizer(max_df=0.90, min_df=5, stop_words='english')
tf = tf_vectorizer.fit_transform(X)
len(tf_vectorizer.get_feature_names())

# LDAのモデル作成と学習
lda = LatentDirichletAllocation(
                            n_components=7,
                            learning_method='online',
                            max_iter=20
                        )
lda.fit(tf)

それでは、学習した7個のトピックについて、それぞれの頻出語をみてみます。


features = tf_vectorizer.get_feature_names()

for tn in range(7):
    print("topic #"+str(tn))
    row = lda.components_[tn]
    words = ', '.join([features[i] for i in row.argsort()[:-20-1:-1]])
    print(words, "\n")

出力は下記の通りです。(乱数の影響で、モデルの学習をやり直すと結果は変わります。)

topic #0
god, people, think, don, know, just, like, does, say, believe, jesus, church, time, way, did, christ, things, good, christian, question

topic #1
25, 10, 11, 12, 14, 16, 15, 17, 20, 13, 18, 19, 55, 30, la, period, 24, 21, pit, 92

topic #2
armenian, armenians, turkish, people, turkey, armenia, turks, greek, genocide, russian, azerbaijan, government, history, muslim, university, soviet, war, 000, ottoman, killed

topic #3
game, don, said, team, just, didn, hockey, like, know, went, year, time, games, think, got, people, going, did, ll, came

topic #4
israel, jews, jewish, israeli, arab, state, people, world, right, public, arabs, rights, human, war, anti, peace, adl, states, country, palestinian

topic #5
medical, health, disease, cancer, patients, use, new, hiv, doctor, season, good, treatment, years, aids, high, drug, number, time, information, vitamin

topic #6
edu, use, file, window, com, server, program, dos, windows, available, motif, using, version, widget, sun, set, display, mit, x11, information

#1があまり意味のない数値を引き受けてくれていますが、
それ以外は、トピックごとに、宗教や国際的な話題、スポーツに医療に、コンピューターなどの単語が分類されています。

最後に、元の各テキストが、カテゴリーごとに妥当なトピック(話題)を持つと判定さているのか可視化してみてみましょう。
どんな可視化方法が一番わかりやすいか色々試したのですが、カラーマップが比較的良さそうでしたので紹介します。
(このほか箱ひげ図などもそこそこ綺麗に特徴が出ましたが。)


topic_data = lda.transform(tf)
fig = plt.figure(figsize=(6, 25))
for i in range(5):
    ax = fig.add_subplot(6, 1, 1+i)
    im = ax.pcolor(topic_data[twenty_news.target == i], vmax=1, vmin=0)
    fig.colorbar(im)
    # 軸の設定
    ax.set_xticks(np.arange(7) + 0.5, minor=False)
    ax.set_xticklabels(np.arange(7))
    ax.set_title(twenty_news.target_names[i])
plt.show()

出力がこちら。

概ね、カテゴリーごとに別のトピックに分類されているのがみて取れます。

投稿日:

scikit-learnのニュースデータをロジスティック回帰で分類

以前書いた、ニュース記事のテキストのサンプルデータを読み込む記事ですが、読み込んだ後何かした記事をまだ書いてなかったのでちょっとやってみようと思います。
参考:20ニュースグループのテキストデータを読み込んでみる
といっても、凝ったモデルを作るのではなく、対した前処理もせずに単純なBoWと ロジスティック回帰だけでどの程度の性能が出るものなか見てみます。

まずはライブラリのインポートと、データの読み込みです。
20種類全部使うとデータが多く時間がかかるのと、miscのカテゴリーも入りどうやってもある程度以上の性能は出ないので、
適当に5カテゴリー選択しました。(categoriesという引数で指定しています。)

また、テキストにヘッダーやフッターが含まれると、ほとんどそこだけの情報で分類できてしまうので、
removeを使って除去しています。
この辺りの仕様はドキュメント参照。
sklearn.datasets.fetch_20newsgroups


# 必要ライブラリのインポート
from sklearn.datasets import fetch_20newsgroups
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from hyperopt import hp,  fmin, tpe, Trials, space_eval

# データの読み込み

remove = ("headers", "footers", "quotes")
categorys = [
        "rec.sport.hockey",
        "soc.religion.christian",
        "sci.med",
        "comp.windows.x",
        "talk.politics.mideast",
    ]

twenty_train = fetch_20newsgroups(
                                subset="train",
                                remove=remove,
                                categories=categorys
                            )
twenty_test = fetch_20newsgroups(
                                subset="test",
                                remove=remove,
                                categories=categorys
                            )

X_train = twenty_train.data
y_train = twenty_train.target
X_test = twenty_test.data
y_test = twenty_test.target

# パラーメータチューニングのため訓練データを2つに分ける
X_train_A, X_train_B, y_train_A, y_train_B = train_test_split(
                                                            X_train,
                                                            y_train,
                                                            test_size=0.2,
                                                            stratify=y_train
                                                        )

見ての通り、X_train, X_test をさらに二つのグループに分けています。
ハイパーパラメーターをチューニングする際に、ここで作ったAグループで学習して、Bグループで評価するようにし、
デフォルトで用意されているテストデータは最後の評価時まで触らずに取っておきます。
本当はクロスバリデーションなどを真面目にやった方がいいのですが、今回はこの方針で。

続いて、ハイパーパラメーターの決定に進みます。
グリッドサーチでもいいのですが、せっかくなのでhyperoptを使ってみます。

spaceは最初はもっと広い範囲で探索していましたが、何度かトライして絞り込みました。
(lr__penalty の l1なども試していたのですが、 Cが大きい時に非常に時間がかかったのと、
今回のデータではl2の方が性能が良かったので探索範囲からも除外。)


# ハイパーパラメーターの探索の準備
space = {
        'cv__min_df': 1 + hp.randint('min_df', 20),
        'cv__max_df': hp.uniform('max_df', 0.5, 0.9),
        'lr__penalty': hp.choice('penalty', ['l2']),
        'lr__C': hp.loguniform('C', -5, 5),
    }


def create_model(args):
    clf = Pipeline(
        [
            ("cv", CountVectorizer(min_df=args["cv__min_df"], max_df=args["cv__max_df"])),
            ("lr", LogisticRegression(C=args["lr__C"], penalty=args["lr__penalty"]))
        ]
    )
    return clf


def objective(args):
    clf = create_model(args)
    clf.fit(X_train_A, y_train_A)
    return - accuracy_score(y_train_B, clf.predict(X_train_B))


trials = Trials()
best = fmin(
            fn=objective,
            space=space,
            algo=tpe.suggest,
            max_evals=100,
            trials=trials,
        )

# 結果の表示
print(best)
print(space_eval(space, best))
# 以下出力
'''
100%|██████████| 100/100 [01:42<00:00,  1.02it/s, best loss: -0.9016949152542373]
{'C': 0.11378962521452059, 'max_df': 0.5208423432316657, 'min_df': 0, 'penalty': 0}
{'cv__max_df': 0.5208423432316657, 'cv__min_df': 1, 'lr__C': 0.11378962521452059, 'lr__penalty': 'l2'}
'''

これでパラーメーターが決まりましたので、モデルを作って評価します。
ここでの評価には最初にとっておいたテストデータを使います。


# 最良のパラメーターでモデルを構築し学習
clf = create_model(space_eval(space, best))
# 訓練データはすべて使う
clf.fit(X_train, y_train)

# 評価
print("正解率:", accuracy_score(y_test, clf.predict(X_test)))
print(classification_report(y_test, clf.predict(X_test), target_names=twenty_test.target_names))

# 以下出力
'''
正解率: 0.8615071283095723
                        precision    recall  f1-score   support

        comp.windows.x       0.90      0.90      0.90       395
      rec.sport.hockey       0.79      0.95      0.86       399
               sci.med       0.87      0.80      0.84       396
soc.religion.christian       0.88      0.82      0.85       398
 talk.politics.mideast       0.87      0.83      0.85       376

           avg / total       0.86      0.86      0.86      1964
'''

正解率約86%。
ほとんど何も工夫していないロジスティック回帰にしてはそこそこの結果だと思います。

投稿日:

scipyの階層的クラスタリングで使える距離関数について

再び階層的クラスタリングの記事の続きです。
参考:scipyで階層的クラスタリング

この記事中では、元のデータの点と点の距離はユークリッド距離を採用しました。
metric='euclidean'と指定しているところがそれです。


# ユークリッド距離とウォード法を使用してクラスタリング
z = linkage(X, metric='euclidean', method='ward')

しかし実際は、これ以外にも多くの距離関数が利用できます。
ドキュメントを読むと、 pdistのページを読めと書いてあるのでそちらをみてみましょう。
scipy.cluster.hierarchy.linkage

metric : str, optional
The distance metric to use. See the distance.pdist function for a list of valid distance metrics.

ちなみに method のほうで指定できる値とその意味は、ちゃんとlinkageのページに載っています。
’single’/’complete’/’average’/’weighted’/’centroid’/’median’/’ward’

scipy.spatial.distance.pdist

metric : str or function, optional
The distance metric to use. The distance function can be ‘braycurtis’, ‘canberra’, ‘chebyshev’, ‘cityblock’, ‘correlation’, ‘cosine’, ‘dice’, ‘euclidean’, ‘hamming’, ‘jaccard’, ‘jensenshannon’, ‘kulsinski’, ‘mahalanobis’, ‘matching’, ‘minkowski’, ‘rogerstanimoto’, ‘russellrao’, ‘seuclidean’, ‘sokalmichener’, ‘sokalsneath’, ‘sqeuclidean’, ‘yule’.

この他、自分で定義した関数も使えるようです。

Y = pdist(X, f)
Computes the distance between all pairs of vectors in X using the user supplied 2-arity function f. For example, Euclidean distance between the vectors could be computed as follows:

dm = pdist(X, lambda u, v: np.sqrt(((u-v)**2).sum()))

それぞれの距離関数はこちらのページから辿っていくと確認しやすいです。
Distance computations (scipy.spatial.distance)

投稿日:

scipyのlinkage関数の結果について

前回の記事の続きです。
参考:scipyで階層的クラスタリング

前回の記事で階層的クラスタリングを実行し可視化するところまで紹介しましたが、
今回は一歩戻ってlinkage関数の戻り値の中身を見てみます。
とりあえず、 linkage matrix をprintして結果を見てみましょう。


from sklearn.datasets import load_iris
from scipy.cluster.hierarchy import linkage
X = load_iris().data[::10, 2:4]
print(X.shape)  # (15, 2)
# ユークリッド距離とウォード法を使用してクラスタリング
z = linkage(X, metric='euclidean', method='ward')
print(z.shape)  # (14, 4)
print(z)
# 以下出力
[[ 2.          3.          0.1         2.        ]
 [ 0.          1.          0.1         2.        ]
 [12.         14.          0.14142136  2.        ]
 [ 4.         16.          0.2081666   3.        ]
 [ 6.          8.          0.31622777  2.        ]
 [ 5.          9.          0.36055513  2.        ]
 [ 7.         11.          0.36055513  2.        ]
 [15.         18.          0.39072582  5.        ]
 [10.         17.          0.43969687  3.        ]
 [13.         23.          0.73598007  4.        ]
 [20.         21.          1.0198039   4.        ]
 [19.         25.          2.00831604  6.        ]
 [24.         26.          3.72312593 10.        ]
 [22.         27.          9.80221064 15.        ]]

前回の記事で可視化したのと同じデータなので、以降の説明は前回の記事中の図と見比べながら読むとわかりやすいと思います。
結果のlinkage matrixは、z.shape の値から分かる通り、14行4列の行列の形をしています。
で、この14という値は、元のデータの個数15個から1減らした値です。
階層的クラスタリングのプロセスの中で、1個ずつグルーピングして集約し、もともと15個あったグループを1つにまとめるのでこうなってます。

そして、列ですが、pythonのインデックスいうと 0列目と1列目はあたらに同じグループに含まれるデータのインデックス、
2列目はそれらの要素orクラスタ間の距離、３列めはそこで新たに作られたクラスタに含まれれるデータの個数を意味します。

具体的には、次のデータは、X[2]とX[3]の距離が0.1で、この二つをまとめて要素が２個のクラスタを作ったことを意味します。
[ 2. 3. 0.1 2. ]
そして、明示はされていませんが、その新しいクラスタには、インデックス15が振られます。(元のデータが0~14の15個なので。)

同様に、次のデータで0と1がまとめられてインデックス16のクラスタが作られます。
[ 0. 1. 0.1 2. ]
で、このインデックス16のクラスタは次のデータで4番目の要素とグルーピングされて、要素数3個のクラスタになります。
[ 4. 16. 0.2081666 3. ]
前回の記事のデンドログラムで確かに0と1でできたクラスタに4が合流しているのが描かれていますね。

このようにして、 linkage matrix の中身を直接読むことができます。

投稿日:

scipyで階層的クラスタリング

今回紹介するのは階層型クラスタリングをscipyで実施する方法です。

階層型クラスタリングの各種アルゴリズム自体は、まだエンジニアを本職にしてたころに知り、その時はこれは面白い手法だと感心していたのですが、
いざデータサイエンティストに転職してからはあまり使ってきませんでした。
クラスタリングする時はk-meansなど、クラスタ数を指定する手法を使うことが多いし、
どれとどれが近いとか言った分析は距離行列眺めたり、次元削減してプロットしたりすることが多かったので。
ただ、他の職種のメンバーに説明するときの樹形図(デンドログラム)のわかりやすさを活用したくなり最近使い始めています

さて、本題に戻ります。

階層型クラスタリングを雑に説明すると、一旦個々のデータを全部別々のクラスタに分類し、
その中から近いものを順番に一つのクラスタにまとめるという操作を繰り返し、最終的に全データを1個のクラスタにまとめる方法です。
この操作を途中で打ち切ることで、任意の個数のクラスタにここのデータを分類することができます。
この際、個々の要素の距離をどう定義するのか、またクラスタと要素、クラスタとクラスタの距離をどのように定義するかによって、
手法が複数存在し、結果も変わります。

この階層型クラスタリングを行う関数が、scipyに用意されています。
Hierarchical clustering (scipy.cluster.hierarchy)

非常に多くの関数がありますが使うのは次の3つです。

scipy.cluster.hierarchy.linkage ・・・ 主役。これが階層型クラスタリングの実装。
scipy.cluster.hierarchy.fcluster ・・・ 各要素について、クラスタリングの結果どのクラスタに属するのかを取得する。
scipy.cluster.hierarchy.dendrogram ・・・ 樹形図(デンドログラム)を描く。

ここまで読んであまりイメージがつかめないと思うので、とりあえずやってみましょう。
データは何でもいいのですが、いつものirisでやります。
(そんなに多くの件数も必要ないので、1/10の件数のデータに絞って使い、4つある特徴量のうち2個だけ使います)


# データ取得。 (15件だけ取得し、特徴量も petal lengthとpetal width に絞る
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt
X = load_iris().data[::10, 2:4]

# データを可視化。
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(1, 1, 1, title="iris (sample)")
plt.scatter(X[:, 0], X[:, 1])
for i, element in enumerate(X):
    plt.text(element[0]+0.02, element[1]+0.02, i)
plt.show()

出力がこちら。これをクラスタリングしていきます。
散布図には番号を振っておきましたが、この番号が結果に出てくる樹形図(デンドログラム)内の番号に対応します。

つぎに階層型クラスタリングを実行して可視化します。


# 階層型クラスタリングに使用する関数インポート
from scipy.cluster.hierarchy import linkage
from scipy.cluster.hierarchy import dendrogram

# ユークリッド距離とウォード法を使用してクラスタリング
z = linkage(X, metric='euclidean', method='ward')

# 結果を可視化
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1, 1, 1, title="樹形図")
dendrogram(z)
plt.show()

結果がこちら。

0,1,2,3,4 が早い段階で一つのクラスタにまとまっていたり、 12と14、6と8が早々にまとまっていたりと納得性のある形になっています。

あとは、もっと扱いやすい形で、何番のデータが何番目のクラスタに所属するのかのリストを作りましょう。
ここで fcluster 関数を使います。
詳しくはドキュメントにありますが、 criterion に ‘maxclust’を指定して、最大クラスタ数で決めたり、
criterion に’distance’ を指定して、距離で閾値を指定したりできます。

ドキュメントを読むよりやってみたほうがわかりやすいと思うのでそれぞれやってみます。
(この辺の話は専用に別記事を用意して取り上げるかも。)


from scipy.cluster.hierarchy import fcluster

# クラスタ数を指定してクラスタリング
clusters = fcluster(z, t=3, criterion='maxclust')
for i, c in enumerate(clusters):
    print(i, c)

# 以下出力
0 1
1 1
2 1
3 1
4 1
5 3
6 3
7 3
8 3
9 3
10 2
11 3
12 2
13 2
14 2

0,1,2,3,4 と 5,6,7,8,9,11 と 10,12,13,14 の3グループにきちんと別れました。

この他、樹形樹に横線を引いてその位置で分けることもできます。
距離=3くらいで分ければ同じ３グループに分かれるのですが、せっかくなので別のところで切りましょう。
距離=1.7を閾値にして、4グループに分かれるのを確認します。


# 距離の閾値を決めてクラスタリング
clusters1 = fcluster(z, 1.7, criterion='distance')
for i, c in enumerate(clusters1):
    print(i, c)

# 以下出力
0 1
1 1
2 1
3 1
4 1
5 4
6 3
7 4
8 3
9 4
10 2
11 4
12 2
13 2
14 2

クラスタ1や、クラスタ2は先ほどと同じですが、6,8 と 5,7,9,11 が別のクラスタに別れました。
これらが分かれる理由は上の樹形図を見ていただければ理解できると思います。

投稿日:

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")

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

投稿日:

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()

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

投稿日:

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

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