標準化レーベンシュタイン距離

以前の記事で、レーベンシュタイン距離を計算できるライブラリを紹介しました。
参考:pythonで編集距離(レーベンシュタイン距離)を求める

どちらかというと、ライブラリよりもアルゴリズム側の説明の続きなのですが、
標準化されたレーベンシュタイン距離(normalized Levenshtein distance)というものも提案されています。
これは、二つの文字列のレーベンシュタイン距離を、文字数が多い方の文字数で割った値として定義されます。
固有名詞の名寄せなどでレーベンシュタイン距離を使う場合、
こちらを使った方がうまく行くことが多いようです(個人的な経験から。)

以前紹介した、 python-Levenshtein
実装されているんじゃないかと期待してしばらく調べていたのですが、どうやらこれは実装されてないようです。
特にLevenshtein.ratio という関数に期待したのですがこれは全然違いました。

ということで自分で実装しましょう。


import Levenshtein


def normalized_distance(text0, text1):
    dist = Levenshtein.distance(text0, text1)
    max_len = max(len(text0), len(text1))
    return dist / max_len

関数名は normalized_levenshtein_distance にしたかったが流石に長すぎるので少し短縮。
ただ、これでも長いですね。(pep8的につらい)

これによって、例えば通常のレーベンシュタイン距離では、
“バニラ” と “アイス” の組み合わせも “チョコレート” と “チョコレートアイス” もどちらも距離は3でしたが、
標準化した距離を使うことで、
前者は距離1、後者は 1/3(=0.333…)と算出されるようになり、前者の方が離れてると見なせます。

使ってみた結果がこちら。


print(Levenshtein.distance("アイス", "バニラ"))  # 3
print(Levenshtein.distance("チョコレートアイス", "チョコレート"))  # 3
print(normalized_distance("アイス", "バニラ"))  # 1.0
print(normalized_distance("チョコレートアイス", "チョコレート"))  # 0.3333333333333333

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

出力がこちら。

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

プライバシーポリシーページを更新しました

このブログのフッターからアクセスできるプライバシーポリシーページを実態に合わせて更新しました。
もともとwordpressがデフォルトで作成した英語のページがあったのですが、
今回日本語で作っています。

手順

  1. wordpressの管理画面の左ペインの設定内にあるプライバシーを選択
  2. プライバシーポリシーページを変更する の 新規ページを作成
  3. 作成されたページを実態に合わせて編集

初期設定のページは英語だったのですが、サイト設定を日本語にしているせいか、新規作成したページは日本語でした。
ただ、内容がかなり充実しすぎていて関係ない記載も多々あったのでそれらは削除し、
コメントに関する記載だけ残しています。

そこに、アクセス解析のため入れているGAのことを記載して完了です。

Google アナリティクス利用規約
7. プライバシー に明記されていルルールにちゃんと対応できてないのを気にしていたのですが、これで一安心です。

お客様はプライバシー ポリシーを公開し、そのプライバシー ポリシーで、お客様がデータ収集のために Cookie を使用していることを必ず通知するものとします。

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%。
ほとんど何も工夫していないロジスティック回帰にしてはそこそこの結果だと思います。

NLTKを使う準備をする

普段、文章の形態素解析にはMeCabを使用しているのですが、
とあるサンプルコードを動かそうとした時に、その中でNLTKが使われており、思ったように動かなかったのでそのメモです。

ちなみに、 Anaconda で環境を作ったので、 nltk自体はインストールされていました。


~$ python --version
Python 3.6.8 :: Anaconda, Inc.
~$ pip freeze | grep nltk
nltk==3.3

サンプルコードを動かした時にデータエラーがこちら


LookupError:
**********************************************************************
  Resource punkt not found.
  Please use the NLTK Downloader to obtain the resource:

  >>> import nltk
  >>> nltk.download('punkt')

ドキュメントを読むと、何かダウンロードをやらないといけないようです。

Installing NLTK Data

こういう、基本的な処理であってもエラーになります。


>>> import nltk
>>> nltk.word_tokenize('hello nltk!')
Traceback (most recent call last):
    (略)
LookupError:
**********************************************************************
  Resource punkt not found.
  Please use the NLTK Downloader to obtain the resource:

  >>> import nltk
  >>> nltk.download('punkt')

  Searched in:
    (略)
**********************************************************************

エラーメッセージを見る限りでは、 ‘punkt’ってのだけで良さそうですが、一気に入れてしまっておきましょう。


~$ python
>>> import nltk
>>> nltk.download()
showing info https://raw.githubusercontent.com/nltk/nltk_data/gh-pages/index.xml

CLI内で完結すると思ったら windowが立ち上がったので少しびっくりしました。

all を選んで Downloadします。
結構時間がかかります。

これで動くようになりました。


>>> import nltk
>>> nltk.word_tokenize("hello nltk!")
['hello', 'nltk', '!']

scipyで距離行列を計算する

前の記事でちらっと pdist関数が登場したので、scipyで距離行列を求める方法を紹介しておこうと思います。

距離行列の説明はwikipediaにあります。
距離行列 – Wikipedia

要するに、N個のデータに対して、(i, j)成分がi番目の要素とj番目の要素の距離になっているN*N正方行列のことです。

これは次の二つの関数を使って計算できます。
scipy.spatial.distance.pdist
scipy.spatial.distance.squareform

numpyで適当に5点取ってやってみましょう。


import numpy as np
from scipy.spatial.distance import pdist
from scipy.spatial.distance import squareform

# 出力する桁数を抑える
np.set_printoptions(precision=3)
# 乱数生成
X = np.random.randint(-5, 6, size=(5, 2))
print(X)
'''
[[-2 -4]
 [-3 -4]
 [ 2 -1]
 [ 4 -2]
 [-1 -2]]
'''

y = pdist(X)
print(y)
'''
[1.    5.    6.325 2.236 5.831 7.28  2.828 2.236 3.162 5.   ]
'''

M = squareform(y)
print(M)
'''
[[0.    1.    5.    6.325 2.236]
 [1.    0.    5.831 7.28  2.828]
 [5.    5.831 0.    2.236 3.162]
 [6.325 7.28  2.236 0.    5.   ]
 [2.236 2.828 3.162 5.    0.   ]]
'''

変数Mに距離行列が入っているのが確認できますね。
(1,2)成分が [-2 -4] と [-3 -4] の距離の1,
(1,3)成分が [-2 -4] と [ 2 -1] の距離の5など、きちんと距離が入っています。

なお、 pdistの戻り値自体は、正方行列の形をしておらず、squareformで正方行列に変形する必要があるので注意です。
pdist の戻り値(コード中では変数y)には、上三角行列部分の値が、1行めから順番に入っています。
配列の長さは$(N*(N-1)/2) = 5*4/2 = 10$ です。

pdist 関数は metric という引数に様々な距離関数の値をとることができ、いろんな種類の距離行列を作ることができます。
(デフォルトはユークリッド距離)

ここまで紹介した、pdistは1つのデータセットに対して距離行列を生成しましたが、
cdistという関数を使うと、2個のデータセットからそれぞれ1個ずつ要素を抽出して距離行列を作ることもできます。

こちらはsquareformを使わなくても初めから行列の形で結果が得られます。


from scipy.spatial.distance import cdist
XA = X[:3]
XB = X[3:]
print(cdist(XA, XB))

# 出力
'''
[[6.325 2.236]
 [7.28  2.828]
 [2.236 3.162]]
'''

最初の例においても、
pdistとsquareform を使うよりも、 cdist(X, X) した方が便利かもしれません。
結果は同じになります。

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クラスタ間の距離、3列めはそこで新たに作られたクラスタに含まれれるデータの個数を意味します。

具体的には、次のデータは、X[2]とX[3]の距離が0.1で、この二つをまとめて要素が2個のクラスタを作ったことを意味します。
[ 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くらいで分ければ同じ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 が別のクラスタに別れました。
これらが分かれる理由は上の樹形図を見ていただければ理解できると思います。

pandasでカテゴリー変数を数値に変換する

pandasの関数を使ってカテゴリー変数を数値に変換する方法のメモです。
one-hot表現とは異なり、[“a”, “b”, “c”] みたいなデータを、 “a”は0,”b”は1,”c”は2みたいな感じで数字に変換します。

scikit-learnのデータセットにおける、正解ラベルのようなデータ型にする、と言ったほうがわかりやすいかもしれません。
Rにはカテゴリカル変数の型があるのでそちらの方がイメージしやすい人もいるかも。

使う関数は、pandas.factorizeです。

この関数に変換したいリストを渡すと、数値に変換したリストと、何番めの数がもともと何の値だったのかを示す配列を返してくれます。

実際にやってみましょう。
numpyでランダムに10個のデータを作って、それを数値化します。


import numpy as np
import pandas as pd

data = np.random.choice(["a", "b", "c"], 10)
print(data)
# ['b' 'c' 'a' 'c' 'c' 'c' 'c' 'a' 'b' 'a']
labels, uniques = pd.factorize(data)
print(labels)
# [0 1 2 1 1 1 1 2 0 2]
print(uniques)
# ['b' 'c' 'a']

うまくいきました。
uniques の値をアルファベット順にしたい時は、
factorize する時に sort=Trueを指定するとできます。


labels, uniques = pd.factorize(data, sort=True)
print(labels)
# [1 2 0 2 2 2 2 0 1 0]
print(uniques)
# ['a' 'b' 'c']

できました。