カテゴリー: 機械学習

前回、scikit-learnのトピックモデル(LDA)における評価指標として、Perplexityを算出する方法を紹介しました。
参考: トピックモデルの評価指標Perplexityの実験

今回はgensim版です。
gensimのLDAモデルには log_perplexity と言うメソッドがあるので、それを使うだけです、って話であれば前の記事とまとめてしまってよかったのですが、
話はそう単純では無いので記事を分けました。

さて、 log_perplexity ってメソッドですが、いかにも perplexity の自然対数を返してくれそうなメソッドです。
perplexity が欲しかったら、 $\exp(log\_perplexity)$ を計算すれば良さそうに見えます。
しかし、 log_perplexity は perplexity の自然対数では無いと言う事実を確認できるのが次の実験です。

前回の記事と同じく、4つのトピックにそれぞれ5単語を含む架空の言語で実験します。

import numpy as np
from gensim.corpora.dictionary import Dictionary
from gensim.models import LdaModel


word_list = [
    ["white", "black", "red", "green", "blue"],
    ["dog", "cat", "fish", "bird", "rabbit"],
    ["apple", "banana", "lemon", "orange", "melon"],
    ["Japan", "America", "China", "England", "France"],
]
corpus_list = [
    np.random.choice(word_list[topic], 100)
    for topic in range(len(word_list)) for i in range(100)
]

# 単語と単語IDを対応させる辞書の作成
dictionary = Dictionary(corpus_list)
# LdaModelが読み込めるBoW形式に変換
corpus = [dictionary.doc2bow(text) for text in corpus_list]

# トピック数4を指定して学習
lda = LdaModel(corpus, num_topics=4, id2word=dictionary)

# log_perplexity を出力
print(lda.log_perplexity(corpus))
# -2.173078593289852

出力が $-2.17\dots$です。
正常に学習できていれば、Perplexityは約5になるはずなので、$\log(5)=1.609\dots$が期待されるのに、符号から違います。

ドキュメントをよく読んでみましょう。
log_perplexity

Calculate and return per-word likelihood bound, using a chunk of documents as evaluation corpus.
Also output the calculated statistics, including the perplexity=2^(-bound), to log at INFO level.

これによると、perplexityは$2^{-bound}$だと言うことになっていて、どうやら、
log_perplexity()が返しているのは、boundに相当するようです。

計算してみましょう。

print(2**(-lda.log_perplexity(corpus)))
# 4.509847333880428

正解は5なので、それらしい結果が出ています。
ですがしかし、Perplexityとしてはこの値は良すぎます。
今回のダミーデータで学習している限りは5単語未満に絞り込めるはずがないのです。

実際、モデルが学習した結果を見てみましょう。

print(lda.show_topics(num_words=6))
"""
[
    (0, '0.100*"bird" + 0.098*"dog" + 0.092*"melon" + 0.092*"cat" + 0.089*"orange" + 0.089*"rabbit"'), 
    (1, '0.104*"red" + 0.104*"green" + 0.102*"white" + 0.098*"blue" + 0.092*"black" + 0.084*"fish"'),
    (2, '0.136*"lemon" + 0.134*"apple" + 0.128*"banana" + 0.117*"orange" + 0.116*"melon" + 0.045*"China"'),
    (3, '0.216*"France" + 0.191*"America" + 0.181*"Japan" + 0.172*"England" + 0.163*"China" + 0.011*"apple"')
]
"""

本来は各トピック上位の5単語を0.2ずつの出現確率と予測できていないといけないので、今 Perplexity を計算しているモデルは
そんなに精度が良くないのです。(ハイパーパラメーターのチューニングを何もしてないので。それはそれで問題ですが、今回の議題からは外します。)

おかしいので、ソースコードを眺めてみたのですが、 2を底とする対数を取ってる様子は無く、普通に自然対数が使われていました。
なので、これはドキュメントの誤りとみた方が良さそうです。(将来的に修正されると思います。)

perplexity=e^(-bound)

と考えると、辻褄があいます。

print(np.exp(-lda.log_perplexity(corpus)))
# 8.785288789149925

トピック数を 1〜4 と動かして算出してみると明らかです。
トピック数が1の時は全く絞り込めていないので元の単語数の約20,
2の時は半分に絞れるので約10,
4の時は、ちゃんと学習できたら正解の5（ただし、デフォルトのハイパーパラメーターだとそこまで成功しないのでもう少し大きい値）
が算出されるはずです。

やってみます。

for i in range(1, 7):
    # トピック数を指定してモデルを学習
    lda = LdaModel(corpus, num_topics=i, id2word=dictionary)

    print(f"トピック数: {i}, Perplexity: {np.exp(-lda.log_perplexity(corpus))}")
"""
トピック数: 1, Perplexity: 20.032145913774283
トピック数: 2, Perplexity: 11.33724134037765
トピック数: 3, Perplexity: 8.921203895821304
トピック数: 4, Perplexity: 7.436279264160588
トピック数: 5, Perplexity: 7.558708610631221
トピック数: 6, Perplexity: 5.892976661122544
"""

大体想定通りの結果は出ましたね。

さて、 log_perplexity は perplexity の対数では無く、
perplexity 対数の符号を反転させたものである、と言うのは認識しておかないと大間違いの元になります。

というのも、perplexityは小さい方が良いとされる指標です。
と考えると、log_perplexityの出力が小さいパラメーターを選んでしまいがちですがそれは誤りであることがわかります。
対数を取ってから符号を反転させているので、大きいものを採用しないといけないのですね。

(この他、必ずしもPerplexityが小さければいいモデルだとは言えない、と言う議論もあるのですが、
今日の記事の範囲は超えますので省略します。)

このブログでトピックモデルの記事を書いたことがあるのですが、
トピック数の決め方について書いてないのに気づいたので評価指標を紹介します。

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

トピックモデルのトピック数を決めるときは、Perplexityもしくは、Coherenceと呼ばれる指標を参考にします。
今回の記事では、Perplexityを紹介します。

と言っても、数学的な定義やその意味についてはいろんな場所で紹介されているので、
この記事では趣向を変えて、架空のデータで実験して理解を深めることを目指します。

まず、 Perplexity の定義は、各単語の出現確率(尤度)の逆数の幾何平均です。
(数式はいろんなサイトに乗っているので省略します。
書籍では、奥村学さんの「トピックモデルによる統計的潜在意味解析」などに載っています。)

この定義だけでは意味がわからないのですが、
「分岐数、または選択肢の数を表している」と説明されることが多いです。

例えば、ある文章があって、単語が一つ隠されていたとします。
このとき、LDAによって、その単語の選択肢が2000まで絞り込めていたら、
そのモデルの Perplexity は 2000です。
単語を絞り込めている方が優れたモデルとされるので、この値は低い方が良いモデルです。

まだわかりにくいので、ここから実験をしていきましょう。
次のような架空の世界があったとします。

– その世界の言葉には4個の話題(トピック)がある。
– 各話題ごとに、単語は5個ある。(つまりその世界に単語は20個しか無い)
– 各文章は一つの話題のみに言及する。

(これらの条件は正確にはトピックモデルではなく、混合ユニグラムモデルですが、わかりやすさのためご容赦ください。)

以上の設定のもとで、ランダムに100単語からなる文章をトピックごとに100個生成します。
コードを見ていただけるとわかりますが、４個の話題は色、動物、果物、国です。(なんでも良いのですが。)

import numpy as np

word_list = [
    ["white", "black", "red", "green", "blue"],
    ["dog", "cat", "fish", "bird", "rabbit"],
    ["apple", "banana", "lemon", "orange", "melon"],
    ["Japan", "America", "China", "England", "France"],
]
corpus = [
    " ".join(np.random.choice(word_list[topic], 100))
    for topic in range(len(word_list)) for i in range(100)
]

さて、あとは以前紹介したコードで、LDAモデルを作って、Perplexityを計算してみましょう。
scikit-learnの場合、ドキュメントにある通り、モデルがperplexityというメソッドを持っています。

トピック数はこの例では4が正解だとわかっているので、4を使います。

本当は、データを訓練データと評価データにわけて、評価データでperplexityを計算する必要があるのですが、
今回は実験なので訓練に使ったデータでそのまま評価します。

from sklearn.feature_extraction.text import CountVectorizer
from sklearn.decomposition import LatentDirichletAllocation

tf_vectorizer = CountVectorizer()
bow = tf_vectorizer.fit_transform(corpus)

# LDAのモデル作成と学習
lda = LatentDirichletAllocation(
    n_components=4,
)
lda.fit(bow)
# perplexityの計算
print(lda.perplexity(bow))
# 5.268629755256359

Perplexity は約 5.27 と、 5に近い値が出ましたね。
このLDAモデルで単語が5個くらいまで絞り込めていることがわかります。

Perplexity がトピック数の決定に使えることをみるために、他のトピック数でも計算してみましょう。

for c_num in range(1, 9):
    lda = LatentDirichletAllocation(
        n_components=c_num,
    )
    lda.fit(bow)
    print(f"トピック数: {c_num}, Perplexity: {lda.perplexity(bow)}")
"""
トピック数: 1, Perplexity: 20.033955224623902
トピック数: 2, Perplexity: 10.330848184515682
トピック数: 3, Perplexity: 7.397066706843117
トピック数: 4, Perplexity: 5.268629755256354
トピック数: 5, Perplexity: 5.305381334487885
トピック数: 6, Perplexity: 5.3074106945229875
トピック数: 7, Perplexity: 5.3206895866734305
トピック数: 8, Perplexity: 5.3529382429024315
"""

トピック数が1個の時は、全く絞り込めていないので、全単語数の20に近い値が出ています。
トピック数が2の場合は、半分に絞れているので約10ですね。
そして、トピック数が4の時に、大体5単語に絞れており、
それ以上トピック数を増やしても大きな改善はありません。
このことから、トピック数は4がベストだろうと判断することができます。

現実世界のデータで試すと、こんなに綺麗にトピック数を決めれたことが無く、
Perplexity の有効性に疑問を持っていたのですが、
理論的にはなかなか良い指標であることが確認できました。

※ この記事を書いたときの僕の環境のgensimのバージョンは 3.8.0 です。

gensimでword2vecなり、トピックモデル(LDA)なり、そのための辞書なりを学習したとき、
saveメソッドで学習したモデルを保存し、同様にloadメソッドで読み込んで使うことができます。
このとき、ファイルに保存するので当然ファイル名を決めないといけないのですが、何形式で保存しているのかよくわからず拡張子に悩んでいました。

以前書いた、gensimでword2vec という記事のサンプルコードでは、ファイル形式がわからないのでとりあえず .model としています。
公式サイトの、Usage examples もそうなってるのですよ。

LDAの方は、
gensim.models.ldamodel.LdaModel.save
の Note を読むと、どうやら、pickleを使ってるように見えます。

ソースコードを確認してみましょう。
https://github.com/RaRe-Technologies/gensim/blob/master/gensim/models/ldamodel.py

def save(~略~):
のところを見ていくと、
salf.state ってのを使って、saveしていますね。

def __init__(~略~):
のところを確認すると、salf.state には、 LdaState というクラスのインスタンスが格納されており、
LdaStateは、
class LdaState(utils.SaveLoad):
とある通り、utils.SaveLoadを継承しています。
どうやらこれが保存と読み込みの本体のようです。

そのソースコードがこちらです。
https://github.com/RaRe-Technologies/gensim/blob/master/gensim/utils.py

このファイル内の、
class SaveLoad:
のところを見ていくと、普通にpickleを使ってファイルに書き出したり読み込んだ理されていることがわかりました。

LDAの次はWord2Vecです。
こちらは話が単純です。

ソースコードを見てみます。
https://github.com/RaRe-Technologies/gensim/blob/master/gensim/models/word2vec.py

class Word2Vec(utils.SaveLoad):

というふうに、モデル自体が、utils.SaveLoadを継承して作られており、先ほどのLDAと同様に保存と読み込みにはpickleが使われています。

Dictionary も同様です。
https://github.com/RaRe-Technologies/gensim/blob/master/gensim/corpora/dictionary.py

pickleは、このブログでも以前記事にしたことがあるようにPythonのオブジェクトを手軽にファイルに書き出せる形式です。
参考: pickleを使ってpythonのオブジェクトをファイルに保存する
ただ、これを使うとなると、いちいちwith openとかいろいろ書かないといけなくてややこしいので、
saveやloadなどわかりやすいメソッド名でラッピングしてもらえているのはありがたいですね。

さて、冒頭に挙げた問題の拡張子名ですが、.pickle あたりを使えば良さそうです。

最近久々にgensimのトピックモデルを使う機会がありました。
そのとき、出現する単語を頻度で絞り込みたったので方法を調べました。

トピックモデルの方法自体は、既に記事を書いてますのでこちらをご参照ください。
参考: gensimでトピックモデル(LDA)をやってみる

さて、gensimのLDAは、学習するコーパスを(単語, 出現回数) というタプルの配列に変換して読み込ませる必要があり、
その形への変換に、gensim.corpora.dictionary.Dictionaryを使います。
この辞書は、何も指定しないと、1回以上出現した単語を全部学習してしまいます。
それを、Scikit-learnのCountVectorizerで、min_dfを指定したときみたいに、n回以上出現した単語のみ、と足切りしたいというのが今回の記事の目的です。

Dictionaryの語彙学習時に指定できる引数の中に、CountVectorizerのmin_dfに相当するものがなかったので、てっきり指定できないのかと思っていたのですが、
じつは、学習した後に、語彙を絞り込む関数である、filter_extremesが用意されていることがわかりました。

使いかたを説明するために、まず適当な単語の羅列でコーパスを作って、辞書を学習しておきます。

import numpy as np
from gensim.corpora.dictionary import Dictionary

#　単語リスト作成
words = [
    "White",
    "Black",
    "Grey",
    "Red",
    "Orange",
    "Yellow",
    "Green",
    "Purple",
    "Blue",
    "Cyan",
    "Magenta",
]

# 再現性のためシードを固定する
np.random.seed(2)
# 単語を適当に選んで文章データを生成
documents = [
    np.random.choice(words, size=np.random.randint(3, 7)).tolist() for _ in range(10)
]

print(documents)
"""
[['Blue', 'Green', 'Grey'],
 ['Blue', 'Purple', 'Grey', 'Black', 'Yellow', 'Magenta'],
 ['Orange', 'Yellow', 'Purple'],
 ['Green', 'Orange', 'Magenta', 'Red', 'Purple', 'Green'],
 ['Black', 'Magenta', 'Red', 'Yellow', 'Blue'],
 ['Green', 'Red', 'Cyan'],
 ['Grey', 'White', 'Orange', 'Grey', 'Orange', 'Magenta'],
 ['Black', 'Purple', 'Blue', 'Grey', 'Magenta', 'Cyan'],
 ['Blue', 'Purple', 'Black', 'Green'],
 ['Magenta', 'Yellow', 'Cyan']]
"""

#　辞書の作成
dictionary = Dictionary(documents)

# 学習した単語リスト
for word, id_ in dictionary.token2id.items():
    print(f"id: {id_}, 単語: {word}, 出現ドキュメント数: {dictionary.dfs[id_]}, 出現回数: {dictionary.cfs[id_]}")

"""
id: 0, 単語:  Blue,    出現ドキュメント数: 5, 出現回数: 5
id: 1, 単語:  Green,   出現ドキュメント数: 4, 出現回数: 5
id: 2, 単語:  Grey,    出現ドキュメント数: 4, 出現回数: 5
id: 3, 単語:  Black,   出現ドキュメント数: 4, 出現回数: 4
id: 4, 単語:  Magenta, 出現ドキュメント数: 6, 出現回数: 6
id: 5, 単語:  Purple,  出現ドキュメント数: 5, 出現回数: 5
id: 6, 単語:  Yellow,  出現ドキュメント数: 4, 出現回数: 4
id: 7, 単語:  Orange,  出現ドキュメント数: 3, 出現回数: 4
id: 8, 単語:  Red,     出現ドキュメント数: 3, 出現回数: 3
id: 9, 単語:  Cyan,    出現ドキュメント数: 3, 出現回数: 3
id: 10, 単語: White,   出現ドキュメント数: 1, 出現回数: 1
"""

これを4個以上の文章に登場した単語だけに絞りこみたいとすると、
filter_extremes(no_below=4)
を実行すれば良いよいうに思えます。
それでやってみたのがこちら。

dictionary.filter_extremes(no_below=4)
for word, id_ in dictionary.token2id.items():
    print(f"単語: {word}, 出現ドキュメント数: {dictionary.dfs[id_]}")
"""
単語: Blue, 出現ドキュメント数: 5
単語: Green, 出現ドキュメント数: 4
単語: Grey, 出現ドキュメント数: 4
単語: Black, 出現ドキュメント数: 4
単語: Purple, 出現ドキュメント数: 5
単語: Yellow, 出現ドキュメント数: 4
"""

Orange/Red/Cyan/White が消えましたね。Orangeは出現回数自体は4でしたが、ドキュメント数が3だったので消えています。
ここで注意なのが、出現ドキュメント数が6だった、Magentaも消えていることです。

これは、filter_extremesのデフォルトの引数が、(no_below=5, no_above=0.5, keep_n=100000, keep_tokens=None) と、
no_above=0.5 も指定されていることに起因します。
つまり、全体の0.5=50%よりも多く出現している単語も一緒に消してしまうわけです。
逆に、no_above だけ指定しても、no_belowは5扱いなので、4文書以下にしか登場しない単語は足切りされます。

この挙動が困る場合は、忘れないように、no_belowとno_aboveを両方指定する必要があります。

#　もう一度辞書の作成
dictionary = Dictionary(documents)
dictionary.filter_extremes(no_below=4, no_above=1)
for word, id_ in dictionary.token2id.items():
    print(f"単語: {word}, 出現ドキュメント数: {dictionary.dfs[id_]}")
"""
単語: Blue, 出現ドキュメント数: 5
単語: Green, 出現ドキュメント数: 4
単語: Grey, 出現ドキュメント数: 4
単語: Black, 出現ドキュメント数: 4
単語: Magenta, 出現ドキュメント数: 6
単語: Purple, 出現ドキュメント数: 5
単語: Yellow, 出現ドキュメント数: 4
"""

出現回数で足切りするのではなく、残す単語数を指定したい場合は、keep_n を使えます。
(これにもデフォルト引数が入ってるので気をつけてください。元の単語数が100000を超えていたら、意図せず動作します)

5単語に絞り込むコードはこうなります。
単語は出現頻度が高い順に選ばれます。
no_below や no_aboveも同時に作用するので、これらの設定次第では、keep_nで指定したよりも少ない単語しか残らないことがあります。

#　もう一度辞書の作成
dictionary = Dictionary(documents)
dictionary.filter_extremes(no_below=1, no_above=1, keep_n=5)
for word, id_ in dictionary.token2id.items():
    print(f"単語: {word}, 出現ドキュメント数: {dictionary.dfs[id_]}")
"""
単語: Blue, 出現ドキュメント数: 5
単語: Green, 出現ドキュメント数: 4
単語: Grey, 出現ドキュメント数: 4
単語: Magenta, 出現ドキュメント数: 6
単語: Purple, 出現ドキュメント数: 5
"""

あとは、あまり使わなさそうですが、 keep_tokens に単語を指定することで、no_belowや、no_aboveに関係なく、
その単語を残すことができます。

#　もう一度辞書の作成
dictionary = Dictionary(documents)
dictionary.filter_extremes(no_below=5, no_above=1, keep_tokens=["White"])
for word, id_ in dictionary.token2id.items():
    print(f"単語: {word}, 出現ドキュメント数: {dictionary.dfs[id_]}")
"""
単語: Blue, 出現ドキュメント数: 5
単語: Magenta, 出現ドキュメント数: 6
単語: Purple, 出現ドキュメント数: 5
単語: White, 出現ドキュメント数: 1
"""

小ネタですが、Dictionaryオブジェクトは、各単語が出現したドキュメント数をdfs, 出現した回数をcfsという変数に保有しています。
filter_extremes を実行すると、dfsの方は単語が絞り込まれた上でidも振り直されるのですが、
cfsは単語が絞り込まれるだけで、idが振り直されません。
(なぜこんな仕様になっているのかは謎です。将来的に修正されるような気がします。)
直前のサンプルコードを動かした時点で、 dfsとcfs の中身を見たものがこちらです。
単語数が4個に減っているのは共通ですが、cfsの方はidが5とか10とか、元のままであることがわかります。

print(dictionary.dfs)
# {0: 5, 2: 5, 1: 6, 3: 1}
print(dictionary.cfs)
# {0: 5, 5: 5, 4: 6, 10: 1}

前回の記事で書いたのがNumPyで多項式回帰する方法だったので、今回はscikit-learnで行う方法を紹介します。
参考: NumPyで多項式回帰

NumPyの方法は、1変数の多項式に特化していたので、変数が1個しかない場合は非常に手軽に使えました。
ただ、実際は $x_0$ の多項式に加えて $x_1$, $x_2$ の変数も使って回帰したいとか、
$x_0, x_1$ を組み合わせた $x_0 * x_1$ みたいな項も入れたいとかいろんなケースがあると思います。
そのような場合は、scikit-learnの利用が検討できます。

と言っても、scikit-learnに多項式関数のモデルが実装されているわけではなく、
実際は前処理だけやってくれるモデルと、通常の線形回帰のモデルを組み合わせて使うことになります。
(このめんどくささが、1変数ならNumPyを推す所以です。)

多項式の特徴量生成には、 sklearn.preprocessing.PolynomialFeatures を使います。

試しに、3変数のデータ4セットに対して、2次までの項を生成してみたコードが次です。

import numpy as np
from sklearn.preprocessing import PolynomialFeatures

X = np.arange(12).reshape(4, 3)
print(X)
"""
[[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]]
"""
poly = PolynomialFeatures(2)  # 次数２を指定
X_poly = poly.fit_transform(X)
print(X_poly)
"""
[[  1.   0.   1.   2.   0.   0.   0.   1.   2.   4.]
 [  1.   3.   4.   5.   9.  12.  15.  16.  20.  25.]
 [  1.   6.   7.   8.  36.  42.  48.  49.  56.  64.]
 [  1.   9.  10.  11.  81.  90.  99. 100. 110. 121.]]
 """

元のデータを$x_0, x_1, x_2$ とすると、
定数$1$,$x_0, x_1, x_2, x_0^2, x_0x_1, x_0x_2, x_1^2, x_1x_2, x_2^2$ のデータが生成されているのがわかります。

何番目のデータがどういう演算で生成された項なのか、という情報は、powers_ という属性に保有されています。

print(poly.powers_)
"""
[[0 0 0]
 [1 0 0]
 [0 1 0]
 [0 0 1]
 [2 0 0]
 [1 1 0]
 [1 0 1]
 [0 2 0]
 [0 1 1]
 [0 0 2]]
"""

定数1の項はいらないな、という時は、 include_biasにFalseを指定して、
PolynomialFeatures(2, include_bias=False)とすれば出てきません。

あとは、この生成されたデータを使って回帰分析を行えば、 scikit-learn を用いた多項式回帰の完成です。

scikit-learnで学習した決定木の学習結果を確認するにはライブラリを使うのが便利ですが、
自分でも直接取得してみたかったので方法を調べてみました。

参考:
dtreevizで決定木の可視化
graphvizで決定木を可視化

とりあえず、 iris を学習しておきます。dtreevizの記事とパラメーターを揃えたので、
この後の結果はそちらと見比べていただくとわかりやすいです。
ただし、最初の分岐が2パターンあって乱数でどちらになるか決まるので、運が悪いと結果が変わります。

from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier

iris = load_iris()
clf = DecisionTreeClassifier(min_samples_split=5)
clf.fit(
    iris.data,
    iris.target
)

"""
DecisionTreeClassifier(ccp_alpha=0.0, class_weight=None, criterion='gini',
                       max_depth=None, max_features=None, max_leaf_nodes=None,
                       min_impurity_decrease=0.0, min_impurity_split=None,
                       min_samples_leaf=1, min_samples_split=5,
                       min_weight_fraction_leaf=0.0, presort='deprecated',
                       random_state=None, splitter='best')
"""

ロジスティック回帰などであれば、係数が coef_に入っているだけなので簡単なのですが、
決定木の場合読み解くのに少し手間がかかります。

その辺りのことは、ドキュメントにも
Understanding the decision tree structureとしてまとめてあるのでこちらも参照しながら読み解いてみました。

必要な情報は clf.tree_の属性としてまとまっているので順番に取り出してみます。

# ノードの数
n_nodes = clf.tree_.node_count
print(n_nodes)
# 13

# 各ノードに振り分けられた学習データの数。
node_values = clf.tree_.value

# 各ノードの左の子ノード。　葉の場合は -1
children_left = clf.tree_.children_left
print(children_left)
# [ 1 -1  3  4  5 -1 -1  8 -1 -1 11 -1 -1]

# 各ノードの右の子ノード。 葉の場合は -1
children_right = clf.tree_.children_right
print(children_right)
# [ 2 -1 10  7  6 -1 -1  9 -1 -1 12 -1 -1]

# 分割に使う特徴量。　葉の場合は-2
feature = clf.tree_.feature
print(feature)
# [ 3 -2  3  2  3 -2 -2  3 -2 -2  2 -2 -2]

# 分割に使う閾値。　葉の場合は-2
threshold = clf.tree_.threshold
print(threshold)
"""
[ 0.80000001 -2.          1.75        4.95000005  1.65000004 -2.
 -2.          1.55000001 -2.         -2.          4.85000014 -2.
 -2.        ]
"""

要するに、各ノードが配列の要素に対応しており、
それぞれ配列に、左の子ノード、右の子ノード、分割に使う特徴量、分割に使う閾値が順番に入っています。

これらの情報を日本語に変化して表示すると次の様になるでしょうか。

for i in range(n_nodes):
    print("\nノード番号:", i)
    if children_left[i] == -1:
        print("    このノードは葉です。")
        print("        予測結果: ")
        for v, t in zip(node_values[i][0], iris.target_names):
            print("            "+t+": ", round(v/sum(node_values[i][0]), 3))
    else:
        print(
            "    "+iris.feature_names[feature[i]],
            "が",
            round(threshold[i], 3),
            "未満の場合、ノード:",
            children_left[i],
            "に進み、それ以外の場合は、",
            children_right[i],
            "に進む。"
        )

出力結果のテキストはこちらです。

ノード番号: 0
    petal width (cm) が 0.8 未満の場合、ノード: 1 に進み、それ以外の場合は、 2 に進む。

ノード番号: 1
    このノードは葉です。
        予測結果: 
            setosa:  1.0
            versicolor:  0.0
            virginica:  0.0

ノード番号: 2
    petal width (cm) が 1.75 未満の場合、ノード: 3 に進み、それ以外の場合は、 10 に進む。

ノード番号: 3
    petal length (cm) が 4.95 未満の場合、ノード: 4 に進み、それ以外の場合は、 7 に進む。

ノード番号: 4
    petal width (cm) が 1.65 未満の場合、ノード: 5 に進み、それ以外の場合は、 6 に進む。

ノード番号: 5
    このノードは葉です。
        予測結果: 
            setosa:  0.0
            versicolor:  1.0
            virginica:  0.0

ノード番号: 6
    このノードは葉です。
        予測結果: 
            setosa:  0.0
            versicolor:  0.0
            virginica:  1.0

ノード番号: 7
    petal width (cm) が 1.55 未満の場合、ノード: 8 に進み、それ以外の場合は、 9 に進む。

ノード番号: 8
    このノードは葉です。
        予測結果: 
            setosa:  0.0
            versicolor:  0.0
            virginica:  1.0

ノード番号: 9
    このノードは葉です。
        予測結果: 
            setosa:  0.0
            versicolor:  0.667
            virginica:  0.333

ノード番号: 10
    petal length (cm) が 4.85 未満の場合、ノード: 11 に進み、それ以外の場合は、 12 に進む。

ノード番号: 11
    このノードは葉です。
        予測結果: 
            setosa:  0.0
            versicolor:  0.333
            virginica:  0.667

ノード番号: 12
    このノードは葉です。
        予測結果: 
            setosa:  0.0
            versicolor:  0.0
            virginica:  1.0

先日可視化した結果とバッチリ対応していますね。

※この記事では dtreevizの version 0.8.2 を使っています。

前回の記事では、dtreeviz を使って学習済みの決定木を可視化しました。
dtreevizではこの他にも、1個か2個の特徴量とラベルの関係を可視化できます。

それが、 ctreeviz_univar と、ctreeviz_bivar です。
扱える特徴量がuniの方が１個、biの方が2個です。

データは必要なので、irisを読み込んでおきます。今回は木は不要です。
(その代わり、max_depsかmin_samples_leafのどちらかの設定が必須です。)

from sklearn.datasets import load_iris
iris = load_iris()

まず1個のほうをやってみます。
特徴量4個しかないので全部出します。

import matplotlib.pyplot as plt
from dtreeviz.trees import ctreeviz_univar

figure = plt.figure(figsize=(13, 7), facecolor="w")
for i in range(4):
    ax = figure.add_subplot(2, 2, i+1)
    ctreeviz_univar(
        ax,
        iris.data[:, i],
        iris.target,
        max_depth=2,
        feature_name=iris.feature_names[i],
        class_names=iris.target_names.tolist(),
        target_name='types'
    )

plt.tight_layout()
plt.show()

どの特徴量が有効なのか、自分的にはこれまでで一番わかりやすいと感じました。

次は２個の方です。特徴量２種類とラベルを渡すと、それらの関係を可視化してくれます。
２個ずつ選んで２つのグラフで可視化してみました。
引数、ですがfeature_name　が feature_names になっており、渡す値も文字列が配列になっているので注意が必要です。

from dtreeviz.trees import ctreeviz_bivar

figure = plt.figure(figsize=(5, 12), facecolor="w")
ax = figure.add_subplot(2, 1, 1)
ctreeviz_bivar(
    ax,
    iris.data[:, :2],
    iris.target,
    max_depth=2,
    feature_names=iris.feature_names[:2],
    class_names=iris.target_names.tolist(),
    target_name='types'
)

ax = figure.add_subplot(2, 1, 2)
ctreeviz_bivar(
    ax,
    iris.data[:, 2:],
    iris.target,
    max_depth=2,
    feature_names=iris.feature_names[2:],
    class_names=iris.target_names.tolist(),
    target_name='types'
)

plt.show()

出力がこちら。

これもわかりやすいですね。

早速、前回の記事でインストールした dtreeviz を使ってみます。

※この記事では dtreevizの version 0.8.2 を使っています。
1.0.0 では一部引数の名前などが違う様です。(X_train が x_dataになるなど。)

とりあえず、データと可視化する木がないと話にならないので、いつものirisで作っておきます。

from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier

iris = load_iris()
clf = DecisionTreeClassifier(min_samples_split=5)
clf.fit(
    iris.data,
    iris.target
)

さて、これで学習したモデル(コード中のclf)を可視化します。
リポジトリのコードを見ながらやってみます。

まず、一番シンプルな可視化は、 dtreeviz.trees.dtreevizにモデルと必要なデータを全部渡すものの様です。
(省略不可能な引数だけ設定して実行しましたが、結構多いですね。)

from dtreeviz.trees import dtreeviz

tree_viz = dtreeviz(
    tree_model=clf,
    X_train=iris.data,
    y_train=iris.target,
    feature_names=iris.feature_names,
    target_name="types",
    class_names=iris. target_names.tolist(),
)
tree_viz

出力がこちら。

graphvizで決定木を可視化 でやったのと比べて、とてもスタイリッシュで解釈しやすいですね。

orientation(デフォルトは’TD’)に’LR’を指定すると、向きを縦から横に変更できます。

tree_viz = dtreeviz(
    tree_model=clf,
    X_train=iris.data,
    y_train=iris.target,
    feature_names=iris.feature_names,
    target_name="types",
    class_names=iris. target_names.tolist(),
    orientation='LR',
)
tree_viz

出力がこちら。

木のサイズによってはこれも選択肢に入りそうですね。

本記事の免責事項:
dtreevizの公式ではpipでのインストールが推奨されているようです。
手順を見ると、condaでgraphviz が入っている場合はそのアンインストールまで明記されています。
そのため、本記事を真似される場合は自己責任でお願いします。
Python環境の破損やその他の動作不良の責任は負いません。
自分自身、将来的にそれらの事象が発生したらpipで入れ直す可能性もあります。

また、この記事でインストールしたのは、version 0.8.2 です。
最新のバージョンでは挙動が異なる可能性があります。

免責事項終わり。

さて、決定木をとても綺麗に可視化してくれるという dtreeviz というライブラリがあるのを聞いて以来、試したいと思っていましたが、
conda(と、conda-forge)のリポジトリには見つからないので後回しにしていたのをやってみることにしました。

個人のMacでは環境構築をcondaに統一しているので、pipはあまり使いたくありません。
しかし免責事項の通り、ドキュメントではpipが推奨されています。

自分の端末ならよかろうということで(職場の端末で試す前の毒見も兼ねて)condaで入れることにしました。
使うのは conda skeleton　です。
このブログのこちらの記事が参考になります。
PyPIのパッケージをcondaでインストールする方法

dot や python-graphviz など、必要ライブラリがすでに入っているのもあり、非常にスムーズにインストールできました。

# skeleton で dtreeviz インストール
$ conda skeleton pypi dtreeviz
$ conda build dtreeviz
$ conda install --use-local dtreeviz

# インストール結果確認
$ conda list dtreeviz
# packages in environment at {HOMEPATH}/.pyenv/versions/anaconda3-2019.10:
#
# Name                    Version                   Build  Channel
dtreeviz                  0.8.2            py37h39e3cac_0    local

# 一次ファイル削除
$ conda build purge

Pythonで
from dtreeviz.trees import dtreeviz
をやってみると無事にインポートできたので、導入は成功した様です。
これからの記事で使い方とか書いていきたいと思います。

前回に引き続き、多言語 Universal Sentence Encoder の話です。
テキストをベクトル化しただけで終わるとつまらないので、これを使って、先日のライブドアニュースコーパスの記事分類をやってみました。
最初、本文でやろうとしたのですが、文ベクトルを得るのに結構時間がかかったので、記事タイトルでカテゴリー分類をやってみます。

すごく適当ですが、512次元のベクトルに変換したデータに対してただのニューラルネットワークで学習してみました。

まずはデータの準備からです。

import pandas as pd
import tensorflow_hub as hub
# import numpy as np
import tensorflow_text
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import OneHotEncoder

# ライブドアニュースコーパス データを読み込む
df = pd.read_csv("./livedoor_news_corpus.csv")
# 訓練データと評価データに分割する
df_train, df_test = train_test_split(df, test_size=0.2, stratify=df.category)
df_train = df_train.copy()
df_test = df_test.copy()
df_train.reset_index(inplace=True, drop=True)
df_test.reset_index(inplace=True, drop=True)

# USEモデルの読み込みと、テキストのベクトル化
url = "https://tfhub.dev/google/universal-sentence-encoder-multilingual/3"
embed = hub.load(url)

X_train = embed(df_train.title).numpy()
X_test = embed(df_test.title).numpy()

# 正解ラベル(記事カテゴリ)を One-Hot 表現に変換

ohe = OneHotEncoder()
ohe.fit(df_train.category.values.reshape(-1, 1))
y_train = ohe.transform(df_train.category.values.reshape(-1, 1)).toarray()
y_test = ohe.transform(df_test.category.values.reshape(-1, 1)).toarray()

あとはモデルを作って学習していきます。

from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Input
from tensorflow.keras.layers import Dense
from tensorflow.keras.layers import Dropout
from tensorflow.keras.callbacks import EarlyStopping
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix

# モデルの作成
model = Sequential()
model.add(Input(shape=(512, )))
model.add(Dropout(0.3))
model.add(Dense(128, activation='tanh'))
model.add(Dropout(0.4))
model.add(Dense(32, activation='tanh'))
model.add(Dropout(0.5))
model.add(Dense(9, activation='softmax'))
print(model.summary())
"""
Model: "sequential"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
dropout (Dropout)            (None, 512)               0         
_________________________________________________________________
dense (Dense)                (None, 128)               65664     
_________________________________________________________________
dropout_1 (Dropout)          (None, 128)               0         
_________________________________________________________________
dense_1 (Dense)              (None, 32)                4128      
_________________________________________________________________
dropout_2 (Dropout)          (None, 32)                0         
_________________________________________________________________
dense_2 (Dense)              (None, 9)                 297       
=================================================================
Total params: 70,089
Trainable params: 70,089
Non-trainable params: 0
_________________________________________________________________
"""

model.compile(
    loss="categorical_crossentropy",
    optimizer="adam",
    metrics=['acc']
)

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

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

"""
Train on 5893 samples, validate on 1474 samples
Epoch 1/100
5893/5893 - 3s - loss: 1.9038 - acc: 0.3655 - val_loss: 1.5009 - val_acc: 0.5611
Epoch 2/100
5893/5893 - 1s - loss: 1.3758 - acc: 0.5564 - val_loss: 1.1058 - val_acc: 0.6771

# -- 中略 --

Epoch 28/100
5893/5893 - 1s - loss: 0.7199 - acc: 0.7611 - val_loss: 0.5913 - val_acc: 0.8012
Epoch 29/100
5893/5893 - 1s - loss: 0.7099 - acc: 0.7663 - val_loss: 0.5932 - val_acc: 0.7985
Epoch 30/100
5893/5893 - 1s - loss: 0.7325 - acc: 0.7597 - val_loss: 0.5935 - val_acc: 0.8005
"""

かなり適当なモデルですが、それでもテストデータで80%くらい正解できたようですね。
classification_reportもみておきましょう。

print(classification_report(model.predict_classes(X_test),y_test.argmax(axis=1),target_names=ohe.categories_[0]))
"""
                precision    recall  f1-score   support

dokujo-tsushin       0.77      0.81      0.79       166
  it-life-hack       0.80      0.82      0.81       169
 kaden-channel       0.80      0.80      0.80       174
livedoor-homme       0.56      0.72      0.63        79
   movie-enter       0.85      0.79      0.82       187
        peachy       0.67      0.69      0.68       166
          smax       0.90      0.89      0.89       175
  sports-watch       0.88      0.89      0.89       177
    topic-news       0.88      0.75      0.81       181

      accuracy                           0.80      1474
     macro avg       0.79      0.80      0.79      1474
  weighted avg       0.81      0.80      0.80      1474
"""

どのカテゴリを、どのカテゴリーに間違えたのかを確認したのが次の表です。

df_test["predict_category"] = model.predict_classes(X_test)
df_test["predict_category"] = df_test["predict_category"].apply(lambda x: ohe.categories_[0][x])

print(pd.crosstab(df_test.category, df_test.predict_category).to_html())

predict_category	dokujo-tsushin	it-life-hack	kaden-channel	livedoor-homme	movie-enter	peachy	smax	sports-watch	topic-news
category
dokujo-tsushin	134	2	2	4	5	21	0	1	5
it-life-hack	1	139	13	4	2	5	9	1	0
kaden-channel	1	12	139	4	1	1	6	0	9
livedoor-homme	6	3	5	57	9	14	1	2	5
movie-enter	0	2	1	1	148	8	1	5	8
peachy	21	1	4	7	14	114	1	2	5
smax	0	8	9	0	0	1	156	0	0
sports-watch	1	1	0	2	2	2	0	158	14
topic-news	2	1	1	0	6	0	1	8	135

独女通信とPeachyとか、ITライフハックと家電チャンネルなど、記事タイトルだけだと間違えても仕方がないような誤判定があるくらいで概ね正しそうです。