PythonでF検定を実装する

最近、とある二つのサンプルの分散が異なることを確認する機会があり、F検定を行う必要がありました。
いつもみたいにSciPyですぐにできるだろうと考えていたのですが、
Statistical functions をみる限りでは、SciPyにF検定は実装されていなさそうでした。

その代わり、バートレット検定とルビーン検定が実装されているのですが。
この二つの検定については別途勉強するとして、とりあえずF検定を行いたかったので自分で実装しました。
幸い、F分布自体はSciPyにあるので簡単です。
t検定のメソッドを参考にし、サンプルを二つ渡したらF統計量とp値を返してくれる関数ように実装しました。

※F検定自体の説明は今回は省略します。
興味のあるかたには、東大出版会の統計学入門(いわゆる赤本)の244ページ、12.2.4 母分散の比の検定 あたりが参考になります。


from scipy.stats import f
import numpy as np


def ftest(a, b):
    # 統計量Fの計算
    v1 = np.var(a, ddof=1)
    v2 = np.var(b, ddof=1)
    n1 = len(a)
    n2 = len(b)
    f_value = v1/v2

    # 帰無仮説が正しい場合にFが従う確率分を生成
    f_frozen = f.freeze(dfn=n1-1, dfd=n2-1)

    # 右側
    p1 = f_frozen.sf(f_value)
    # 左側
    p2 = f_frozen.cdf(f_value)
    # 小さい方の2倍がp値
    p_value = min(p1, p2) * 2

    # 統計量Fとp値を返す
    return f_value, p_value

きちんと実装できているかどうか確認するため、統計学入門の練習問題を一つ解いておきます。
252ページの練習問題、 12.2 の iii) が分散の検定なのでやってみます。

次のコードのdata_1とdata_2の分散が等しいというのが帰無仮説で、等しくないというのが対立仮説です。


# 問題文のデータを入力
data_1 = np.array([15.4, 18.3, 16.5, 17.4, 18.9, 17.2, 15.0, 15.7, 17.9, 16.5])
data_2 = np.array([14.2, 15.9, 16.0, 14.0, 17.0, 13.8, 15.2, 14.5, 15.0, 14.4])

# F統計量とp値を計算
f_value, p_value = ftest(data_1, data_2)
print(f"F統計量: {f_value:1.3f}")
# F統計量: 1.564
print(f"p値: {p_value:1.3f}")
# p値: 0.516

F統計量は正解と一致しましたし、帰無仮説が棄却できないのも確認できました。

Jupyter Notebookでインタラクティブに関数を実行する

Juputer Notebookで関数を実行するとき、結果を見ながら何度も引数を変更して再実行したくなることはないでしょうか。
毎回コード中の引数を書き換えて実行しても良いのですが、それをGUIで実行できると便利です。

その便利な機能が、ウィジェットとして用意されているのでそれを紹介します。
使うのは ipywidgets です。
ドキュメントはこちら: Jupyter Widgets
実は細かい設定でGUIをいろいろ作れるのですが、今回はとりあえず最低限度の機能でGUIで操作しながらグラフを描写する関数を実行します。

実行するサンプル関数としてリサージュ曲線を描く関数を用意しました。


import numpy as np
import matplotlib.pyplot as plt


def lissajous_curve(a, b, delta, color="b", title="リサージュ曲線"):

    t = np.linspace(0, 2, 601) * np.pi
    x = np.cos(a * t)
    y = np.sin(b * t + delta)

    fig = plt.figure(figsize=(8, 8))
    ax = fig.add_subplot(111, aspect="equal", title=title)
    ax.plot(x, y, color=color, )
    plt.show()

一応説明しておくと、これは変数$t$を媒介変数とし、
$x=A\cos(at), y=B\sin(bt+\delta)$で表される点をプロットした曲線です。
($A,B$はただ図形の縮尺が変わるだけなので今回のコードは$1$で固定しました。)

さて、通常は$a, b, \delta$や色とタイトルに値を指定して実行することで目当てのグラフを描写しますが、
これらの引数をGUIで指定できる様にします。

方法は簡単で、ipywidgets.interact に先ほどの関数と、指定できる引数の範囲を渡すだけです。
引数の範囲は次の様に指定できます。
– (開始, 終了, ステップ) のタプル(ステップを省略したら1)
– 指定できる値のリスト
– テキスト
– 表示用文字列: 渡す値の辞書


from ipywidgets import interact


interact(
    lissajous_curve,
    a=(1, 10, 1),
    b=(1, 10, 1),
    delta=(0, 2*np.pi, 0.01),
    title="リサージュ曲線",
    # ただの配列で指定することも可能。
    # color=["b", "k", "r", "g", "y", "m", "c"]
    # 表示用文字列: 渡す値の辞書も使える
    color={
        "Blue": "b",
        "Green": "g",
        "Red": "r",
        "Black": "k",
    }
)

実行したのがこちらです。
Blogなので画像キャプチャになってしまっていますが、実際はスライドバーやドロップダウンで操作し、グラフを書き直すことができます。

NumPyのブロードキャストで変換できる型

NumPyを普段使いしてると便利な機能に、ブロードキャストがあります。
これは配列のサイズが違うもの通しを演算するときに、いい感じに小さい方を拡張して演算してくれるものです。

例えば、配列とスカラーの和や、行列とベクトルの和を次の様に計算してくれます。


a = np.array(range(4))
print(a)
# [0 1 2 3]

# 7 を [7, 7, 7, 7] として扱って足してくれる
b = a + 7
print(b)
# [ 7  8  9 10]

c = np.array(range(6)).reshape(2, 3)
print(c)
"""
[[0 1 2]
 [3 4 5]]
"""

d = np.array([5, 5, 5])

# [5, 5, 5] を [[5, 5, 5], [5, 5, 5]] として扱って足してくれる
e = c+d
print(e)
"""
[[ 5  6  7]
 [ 8  9 10]]
 """

本当にいい感じにやってくれるのであまり意識せずに使っていましたが、仕様を正確に把握しておきたかったので改めてドキュメントを読みました。
と言うのも、いつでも動くわけではありませんし、正方行列とベクトルの和の様にどちらの軸にブロードキャストされるか迷うことなどあるからです。


# 動かない例
a = np.array(range(6)).reshape(3, 2)
b = np.array(range(3))
a + b
# ValueError: operands could not be broadcast together with shapes (3,2) (3,)

# 行か列か引き伸ばされる方向がわかりにくい例
c = np.zeros(shape=(3, 3))
d = np.array(range(3))
print(c+d)
"""
[[0. 1. 2.]
 [0. 1. 2.]
 [0. 1. 2.]]
"""

ドキュメントはこちらです。
参考: Broadcasting

長いのですが、基本的なルールは General Broadcasting Rules に書かれてる次の法則だけです。
配列の次元数は後ろから順番に前に向かって比較されます。(長さが違う場合は、短い方に1が追加されて揃えられます。)
そして、それらの値が等しいか、もしくは一方が1であればブロードキャストされます。

ブロードキャストされるときは、値が1だった(もしくは無かった)次元の向きにデータがコピーされて拡張されます。

先ほどのエラーが起きた例で言えば、 (3, 2)次元と (3) 次元の 2と3が比較されてこれが等しくないからエラーになったわけですね。
その次の例についてはまず、
[0, 1, 2] (shapeは(3,))に、次元が一個追加されて、
[[0, 1, 2]] (shapeは(1, 3)) に変換され、それが各行にコピーされたので上の例の様な結果になっています。

先述のシンプルなルールさえ満たしていれば、次の例の様な少々無茶でイメージしにくい配列同士でもブロードキャストされます。


a = np.array(range(21)).reshape(1,1,3,1,7)
b = np.array(range(10)).reshape(2,1,5,1)

print(a.shape)
# (1, 1, 3, 1, 7)
print(b.shape)
# (2, 1, 5, 1)

c = a+b
print(c.shape)
# (1, 2, 3, 5, 7)

(1, 1, 3, 1, 7) と (2, 1, 5, 1) では長さが違うので、後者の方の先頭に1が挿入され
以下の二つになるわけですね。
(1, 1, 3, 1, 7)
(1, 2, 1, 5, 1)
これを順番に見ていくと、ブロードキャストのルールをみたいしているので足りない向きについてはデータがコピーされ、
和がとれているわけです。

pandasのDataFrameのappendは遅い

この記事のタイトルは悩んだのですが、一番伝えたかった内容がそれなので、上記の様になりました。
他のタイトル候補は次の様になります。

– 配列を要素に持つDataFrameの縦横変換
– 高速にDataFrameを作成する方法

要するに、1件ごとに発生するデータを毎回appendしてデータフレームを構成する処理は、
コーデイングの観点ではわかりやすいですが、件数によっては非常に時間がかかります。

僕の場合ですが、DataFrameの要素に配列(もしくはカンマ区切りなどの文字列)が入っていた時に、
それを要素ごとに縦横変換する時によく遭遇します。

例えば以下の感じのDataFrameを


       key                                         value
0    key_0          [value_0, value_1, value_2, value_3]
1    key_1                            [value_4, value_5]
2    key_2                   [value_6, value_7, value_8]
3    key_3       [value_9, value_10, value_11, value_12]
4    key_4                          [value_13, value_14]
# (以下略)

次の様なDataFrameに変換したい場合です。


     key    value
0  key_0  value_0
1  key_0  value_1
2  key_0  value_2
3  key_0  value_3
4  key_1  value_4
5  key_1  value_5
6  key_2  value_6
7  key_2  value_7
8  key_2  value_8
# (以下略)

データ量が少ない場合は、1件ずつ appendしても問題ありません。
次の様に、iterrows で1行1行取り出して、出力先にpivot_dfに追加してくと、素直で理解しやすい処理になります。


pivot_df = pd.DataFrame(columns=["key", "value"])

for _, row in df.iterrows():
    value_list = row["value"]
    for v in value_list:
        pivot_df = pivot_df.append(
            pd.DataFrame(
                {
                    "key": [row.key],
                    "value": [v],
                }
            )
        )

ただし、この処理は対象のデータ量が2倍になれば2倍以上の時間がかかる様になります。
そこで、大量データ(自分の場合は数十万件以上くらい。)を処理する時は別の書き方を使っていました。

つい最近まで、自分の環境が悪いのか、謎のバグが潜んでいるのだと思っていたのですが、
ドキュメントにも行の追加を繰り返すと計算負荷が高くなるからリストに一回入れろと書いてあるのを見つけました。
どうやら仕様だった様です。

引用

Iteratively appending rows to a DataFrame can be more computationally intensive than a single concatenate. A better solution is to append those rows to a list and then concatenate the list with the original DataFrame all at once.

その様なわけで、DataFrame中の配列の縦横変換は自分は以下の様に、一度配列を作って配列に追加を続け、
出来上がった配列をDataFrameに変換して完成としています。


keys = []
values = []

for _, row in df.iterrows():
    keys += [row.key] * len(row.value)
    values += row.value

pivot_df = pd.DataFrame(
    {
        "key": keys,
        "value": values,
    }
)

jupyter notebookの拡張機能を導入する

職場のMacでは以前から導入していたのですが、私物のMacにも入れたのでメモしておきます。

そのままでも十分使いやすい jupyter notebookですが、
jupyter_contrib_nbextensions というライブラリを入れるといろいろな拡張機能が使える様になり、一層便利になります。
自動で pep8 を満たす様に整形してくれたり、入力セルを隠せたり、補完機能が強化されたり、
テーブルがきれいに表示されるなどの設定ができます。

インストールはこちらのページにもある通り、
pipでも condaでも入れることができます。 自分はnotebook自体がcondaで導入したものなので、
condaで入れています。
これを jupyter notebookが起動してない時に実行し、その後、notebookを起動します。


# pipで入れる場合のコマンド
# $ pip install jupyter_contrib_nbextensions
# condaで入れる場合のコマンド
$ conda install -c conda-forge jupyter_contrib_nbextensions

すると、 jupyterの画面に、 Nbextensions というタブが生成されます。

チェックボックスのチェックを外すと設定を変えられるようになり、
あとは好みの設定を選んで入れていくだけです。

実際に notebookを触りながら色々試すと楽しいと思います。

(補足) サイトによっては、利用するため(タブを表示する)にはアクティベーションの操作が必要といったことが書いてありますが、
僕の環境ではインストールしたらタブが勝手に出てきました。
バージョンによって違うのかもしれません。

Pythonで特異値分解する方法(SciPy利用)

最近使う機会があったので、特異値分解について紹介します。
まず、$A$をランク$r$の$m$行$n$列の実行列とします。
(複素行列も考慮したい場合は、この後出てくる転置行列を随伴行列に読み替え、直行行列をユニタリ行列で読み替えてください。)

このとき、$U$は$m\times m$の直行行列、$V$は$n\times n$の直行行列、Sは非対角成分が0の$m \times n$行列を用いて、
$$A = USV^{\top}$$
と分解することを特異値分解と言います。
$S$の対角成分のうち$0$でないもの(特異値)の個数が$A$のランク$r$になります。

さて、この特異値分解をPythonで実行する方法は複数あり、numpy, SciPy, scikit-learnにそれぞれ実装があります。
参考:
numpy.linalg.svd
scipy.linalg.svd
sklearn.decomposition.TruncatedSVD

これらの使い分けですが、機械学習のパイプラインに組み込んだり、可視化が目的の時など次元削減のために利用するのであればscikit-learnがおすすめです。
それ以外の場合は、SciPyが良いでしょう。numpyより設定できるオプションが多く、また、いくつか補助的な便利関数が用意されています。

とりあえず乱数で生成した行列を、SciPyで特異値分解してみます。


import numpy as np
import scipy

m = 5
n = 3

np.random.seed(110)  # 乱数固定
A = np.random.randint(-10, 10, size=(m, n))
print(A)
"""
[[-10  -7   5]
 [  5 -10   6]
 [  6  -4  -8]
 [ -2  -1  -5]
 [-10   7  -9]]
"""

U, S_diags, V_t = scipy.linalg.svd(A)

# U
print(U)
"""
U:
[[ 0.10801327  0.81765612 -0.427367   -0.18878079 -0.31857632]
 [ 0.62638958  0.06527027 -0.28451679  0.07142207  0.71925307]
 [ 0.04632544 -0.55033827 -0.69658511 -0.39810108 -0.226421  ]
 [-0.16944935 -0.04718317 -0.43857367  0.87461141 -0.1084836 ]
 [-0.75173806  0.14859275 -0.24254891 -0.18929111  0.56404698]]
"""

# Sの対角成分:
print(S_diags)
# [19.70238709 15.08068891  9.76671719]

# Vの転置行列:"
print(V_t)
"""
[[ 0.51699558 -0.6241887   0.58575083]
 [-0.83177902 -0.20473932  0.51597042]
 [ 0.20213668  0.75396968  0.62503639]]
"""

変数名で表現してみましたが、 $U$以外の二つは少し癖のある形で返ってきます。
$S$は対角成分しか戻ってきませんし、$V$は転置された状態で戻されます。
再度掛け算して結果を角にする時などは注意が必要です。

さて、$S$は対角成分だけ返ってきましたが、これを$m\times n$行列に変換するには、専用の関数が用意されています。
scipy.linalg.diagsvd


S = scipy.linalg.diagsvd(S_diags, m, n)
print(S)
"""
[[19.70238709  0.          0.        ]
 [ 0.         15.08068891  0.        ]
 [ 0.          0.          9.76671719]
 [ 0.          0.          0.        ]
 [ 0.          0.          0.        ]]
"""

あとは、 $USV^{\top}$計算して、元の$A$になることをみておきましょう。


print(U@S@V_t)
"""
[[-10.  -7.   5.]
 [  5. -10.   6.]
 [  6.  -4.  -8.]
 [ -2.  -1.  -5.]
 [-10.   7.  -9.]]
"""

元の値が復元できましたね。

Pythonのlistをsortするときに指定するkey引数について

先日とある目的で、配列を要素に持つ配列を、各配列の3番目の値を使ってソートするコードが必要になりました。
最初はそのキーになる要素だけ取り出して、 argsortしてどうにか並べ替えようとしたのですが、実はkeyって引数を使えば綺麗に実装できることがわかりました。

ドキュメントは組み込み型のlistのsortメソッドの部分が該当します。

key って引数に 単一の引数を取る関数を渡せば、listの各要素にその関数が適用され、戻り値によってソートされる様です。
なので、配列の配列を3番目の要素でソートするには、lambda 式あたりでそう言う関数を作ってあげれば実現できます。
(インデックスは0始まりなので、実際にはインデックス2を取得します。)


# ソートされるデータ作成
list_0 = [
        [9, 8, 5],
        [0, 8, 3],
        [1, 6, 5],
        [9, 0, 0],
        [4, 9, 3],
        [1, 4, 8],
        [4, 0, 6],
        [0, 3, 5],
        [1, 3, 1],
        [5, 2, 7],
    ]

# 3番目(index 2)の要素でソート
list_0.sort(key=lambda x: x[2])
print(list_0)
# [[9, 0, 0], [1, 3, 1], [0, 8, 3], [4, 9, 3], [9, 8, 5], [1, 6, 5], [0, 3, 5], [4, 0, 6], [5, 2, 7], [1, 4, 8]]

これを応用すれば、辞書の配列を特定のkeyに対応する値でソートするといったことも簡単にできます。
(lambda式に辞書オブジェクトのgetメソッド渡すだけですね。)

もちろん、lambda 式以外にも事前に定義した関数を渡すこともできます。
次の例は、トランプのカード (2〜9とTJQKAで表したランクと、cdhsで表したスートの2文字で表したもの)を強い順にソートするものです。


# カードの強さを返す関数
def card_rank(card):
    rank_dict = {
            "T": 10,
            "J": 11,
            "Q": 12,
            "K": 13,
            "A": 14,
        }
    suit_dict = {
            "s": 3,
            "h": 2,
            "d": 1,
            "c": 0,
        }
    return int(rank_dict.get(card[0], card[0]))*4 + suit_dict[card[1]]


card_list = ["Ks", "7d", "Ah",  "3c", "2s", "Td", "Kc"]

# ソート前の並び
print(card_list)
# ['Ks', '7d', 'Ah', '3c', '2s', 'Td', 'Kc']

# 強い順にソート
card_list.sort(key=card_rank, reverse=True)
print(card_list)
# ['Ah', 'Ks', 'Kc', 'Td', '7d', '3c', '2s']

reverse = True を指定しているのはソートを降順にするためです。

このほかにも str.lowerを使って、アルファベットを小文字に統一してソートするなど、
使い方は色々ありそうです。

matplotlibでグラフのx軸とy軸のメモリの間隔(アスペクト)を揃える

matplotlibを使ってる方はご存知だと思いますが、matplotlibはグラフを綺麗に描写するためにx軸とy軸でそれぞれメモリの間隔をいい感じに調整してくれます。
この縦横の比率をアスペクト比と言うそうです。(Wikipedia: アスペクト比)

ほとんどの場合、自動的に調整してくれるのでただありがたいのですが、描写するものによってはこの比率を揃えたいことがあります。
そうしないと円を書きたかったのに楕円になったりします。
実はこのブログの過去の記事のサンプルコードの中で利用したことがあるのですが、
この間必要になったときに自分でこのブログ内から探せなかったので今回独立した記事にしました。

文章で説明していると分かりにくいので、とりあえず縦横の比率がそろってない例を出します。
シンプルにいつものアヤメのデータから2次元分だけ拝借して、散布図を書きました。


import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

# データ取得。 2次元分だけ利用
X = load_iris().data[:, 2:]
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(1, 1, 1, title="アスペクト未指定")
ax.scatter(X[:, 0], X[:, 1])
plt.show()

見慣れた図ですが x軸における幅1 と y軸における幅1が全然違いますね。

この縦横の比率を揃えるには、 add_subplot で Axes オブジェクトを作るときに aspectで指定するか、
Axes オブジェクトの set_aspect 関数で指定します。
指定できる値は “auto”(これがデフォルト), “equal”(比率を揃える), 数値(指定した比率になる) の3パターンです。

個人的には add_subplot した時点でしていておく方が好きです。ただ、行が長くなるので set_aspect 使う方が pep8を守りやすいとも思ってます。

実際にやってみるとこうなります。


fig = plt.figure(facecolor="w")
ax = fig.add_subplot(1, 1, 1, aspect="equal", title="アスペクト equal")
ax.scatter(X[:, 0], X[:, 1])
plt.show()

比率が揃いましたね。
ただ今回の例だと、 “auto”の方が見やすいのでその点ちょっと失敗したと思いました。

set_aspect を使う時は、 anchor という引数を同時に渡すこともできます。
これは aspect の指定によって、グラフが実際より小さくなったときに、元の領域のどの位置に表示するかをしているものです。
‘C’が中心なのはいいとして、 ‘N’とか’SW’とかやけに分かりにくい値で指定しないといけないなと感じていたのですが、どうやらこれは東西南北をEWSNのアルファベットで表したもののようです。

自分は必要になったことがないのですが、 set_anchor と言う関数のドキュメントに説明がありますのでこちらもご参照ください。

WordCloudの文字の色を明示的に指定する

相変わらずWordCloudの話です。(今回くらいで一旦止めます。)

今回は文字の色を個別に指定します。
前回の記事のコードの抜粋が以下ですが、
colormap 引数で どのような色を使うのかはざっくりと指定できます。(実際の色は単語ごとにランダムに決まる)
これを、「この単語は何色」って具体的に指定しようというのが今回の記事です。


# TF-IDFで Word Cloud作成
wc_1 = WordCloud(
        font_path="/Library/Fonts/ipaexg.ttf",
        width=600,
        height=300,
        prefer_horizontal=1,
        background_color='white',
        include_numbers=True,
        colormap='tab20',
    ).generate_from_frequencies(tfidf_dict)

さて、方法ですが、WordCloudのインスタンスを作る時に、color_func って引数で、色を返す関数を渡すことで実装できます。
ドキュメントの該当ページから引用すると、
次のように word とか font_size とかを受け取る関数を用意しておけばいいようです。

color_funccallable, default=None
Callable with parameters word, font_size, position, orientation, font_path, random_state that returns a PIL color for each word. Overwrites “colormap”.

Using custom colors というページに例もあるので、参考にしながらやってみましょう。

さて、どうやって色をつけるかですが、今回は「品詞」で色を塗り分けることにしてみました。ただ、やってみたら名詞のシェアが高すぎてほぼ全体が同じ色になったので、名詞だけは品詞細分類1までみてわけます。

コードですが、前回の記事で準備したデータを使うので、 形態素分析済みで、TF-IDFも計算できているデータはすでにあるものとします。

まず、単語から品詞を返す関数を作ります。
注: 品詞は本当は文脈にも依存するので、形態素解析したときに取得して保存しておくべきものです。ただ、今回そこは本質でないので、一単語渡したらもう一回MeCabにかけて品詞を返す関数を作りました。
もともと1単語を渡す想定ですが、それがさらに複数の単語にわけられちゃったら1個目の品詞を返します。あまりいい関数ではないですね。
あと、処理の効率化のために関数自体はメモ化しておきます。


from functools import lru_cache


@lru_cache(maxsize=None)
def get_pos(word):
    parsed_lines = tagger.parse(word).split("\n")[:-2]
    features = [l.split('\t')[1] for l in parsed_lines]
    pos = [f.split(',')[0] for f in features]
    pos1 = [f.split(',')[1] for f in features]

    # 名詞の場合は、 品詞細分類1まで返す
    if pos[0] == "名詞":
        return f"{pos[0]}-{pos1[0]}"

    # 名詞以外の場合は 品詞のみ返す
    else:
        return pos[0]

(tagger などは前回の記事のコードの中でインスタンス化しているのでこれだけ実行しても動かないので注意してください)

さて、単語を品詞に変換する関数が得られたので、これを使って、単語に対して品詞に応じた色を戻す関数を作ります。


import matplotlib.cm as cm
import matplotlib.colors as mcolors

# 品詞ごとに整数値を返す辞書を作る
pos_color_index_dict = {}
# カラーマップ指定
cmap = cm.get_cmap("tab20")


# これが単語ごとに色を戻す関数
def pos_color_func(word, font_size, position, orientation, random_state=None,
                   **kwargs):

    # 品詞取得
    pos = get_pos(word)

    # 初登場の品詞の場合は辞書に追加
    if pos not in pos_color_index_dict:
        pos_color_index_dict[pos] = len(pos_color_index_dict)

    color_index = pos_color_index_dict[pos]

    # カラーマーップでrgbに変換
    rgb = cmap(color_index)
    return mcolors.rgb2hex(rgb)

**kwargs が吸収してくれるので、実は font_size とか position とか関数中で使わない変数は引数にも準備しなくていいのですが、
いちおうドキュメントの例を参考に近い形で書いてみました。

これを使って、ワードクラウドを作ります。


# TF-IDFで Word Cloud作成
wc = WordCloud(
        font_path="/Library/Fonts/ipaexg.ttf",
        width=600,
        height=300,
        color_func=pos_color_func,
        prefer_horizontal=1,
        background_color='white',
        include_numbers=True,
    ).generate_from_frequencies(tfidf_dict)
wc.to_image()

変わったのは color_func にさっき作った関数を渡しているのと、 colormap の指定がなくなりました。(指定しても上書きされるので無視されます)
そして出力がこちら。

同じ品詞のものが同色に塗られていますね。

どの品詞が何色なのか、凡例?というか色見本も作ってみました。
たまたま対象のテキストで登場した品詞だけ現れます。


fig = plt.figure(figsize=(6, 8), facecolor="w")
ax = fig.add_subplot(111)
ax.barh(
    range(len(pos_color_index_dict)),
    [5] * len(pos_color_index_dict),
    color=[mcolors.rgb2hex(cmap(i)) for i in range(len(pos_color_index_dict))]
)

ax.set_xticks([])
ax.set_yticks(range(len(pos_color_index_dict)))
ax.set_yticklabels(list(pos_color_index_dict.keys()))
plt.show()

単語の頻度データからWord Cloudを作成する方法

今回も Word Cloud の話です。
前回は見た目の設定を変える話でしたが、今回は読み込ませるデータの話になります。

さて、Word Cloudを作るとき、
generate (もしくは generate_from_text) 関数に、 テキストを渡し、その中の出現回数でサイズを決めました。

しかし実際には、出現回数ではなくもっと別の割合でサイズを調整したいことがあります。
例えば、TF-IDFで重みをつけたい場合とか、トピックモデルのトピック別出現確率のようなもともと割合で与えられたデータを使いたい場合などです。

このような時は “単語: 頻度” の辞書(dict)を作成し、
generate_from_frequencies (もしくは fit_words) に渡すと実行できます。
ドキュメントの API Reference には詳しい説明がないので、ギャラリーの Using frequency をみる方がおすすめです。

今回はサンプルとして、ライブドアニュースコーパスから適当に1記事選んで、通常の単語の出現回数(今までと同じ方法)と、
TF-IDFで重み付けした方法(generate_from_frequenciesを使う)でそれぞれ WordCloudを作ってみました。

今まで、単語を名詞動詞形容詞に絞ったり平仮名だけの単語は外したりしていましたが、
今回は差が分かりやすくなるようにするためにSTOPWORDなしで全単語を含めています。(コード中で該当処理をコメントアウトしました。)
「てにをは」系の頻出語が小さくなってSTOPWORD無しでも分かりやすくなる効果が出ているのが感じられると思います。


import re
import MeCab
import pandas as pd
import unicodedata
import matplotlib.pyplot as plt
from wordcloud import WordCloud
from sklearn.feature_extraction.text import TfidfVectorizer

# 分かち書きの中で使うオブジェクト生成
tagger = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd")
# ひらがなのみの文字列にマッチする正規表現
kana_re = re.compile("^[ぁ-ゖ]+$")


def mecab_tokenizer(text):
    # ユニコード正規化
    text = unicodedata.normalize("NFKC", text)
    # 分かち書き
    parsed_lines = tagger.parse(text).split("\n")[:-2]
    surfaces = [l.split('\t')[0] for l in parsed_lines]
    features = [l.split('\t')[1] for l in parsed_lines]
    # 原型を取得
    bases = [f.split(',')[6] for f in features]
    # 品詞を取得
    # pos = [f.split(',')[0] for f in features]
    # 各単語を原型に変換する
    token_list = [b if b != '*' else s for s, b in zip(surfaces, bases)]
    # 名詞,動詞,形容詞のみに絞り込み
    # target_pos = ["名詞", "動詞", "形容詞"]
    # token_list = [t for t, p in zip(token_list, pos) if (p in target_pos)]
    # ひらがなのみの単語を除く
    # token_list = [t for t in token_list if not kana_re.match(t)]
    # アルファベットを小文字に統一
    token_list = [t.lower() for t in token_list]
    # 半角スペースを挟んで結合する。
    result = " ".join(token_list)
    # 念のためもう一度ユニコード正規化
    result = unicodedata.normalize("NFKC", result)

    return result


# データ読み込み
df = pd.read_csv("./livedoor_news_corpus.csv")
# 形態素解析(全データ)
df["token"] = df.text.apply(mecab_tokenizer)

# サンプルとして用いるテキストを一つ選び、形態素解析する
text = df.iloc[1010].text
tokenized_text = mecab_tokenizer(text)

# tfidfモデルの作成と学習
tfidf_model = TfidfVectorizer(token_pattern='(?u)\\b\\w+\\b', norm=None)
tfidf_model.fit(df.token)

# 対象テキストをtf-idfデータに変換
tfidf_vec = tfidf_model.transform([tokenized_text]).toarray()[0]
# 単語: tf-idfの辞書にする。
tfidf_dict = dict(zip(tfidf_model.get_feature_names(), tfidf_vec))
# 値が正のkeyだけ残す
tfidf_dict = {k: v for k, v in tfidf_dict.items() if v > 0}

# 単語の出現頻度でWord Cloud作成
wc_0 = WordCloud(
        font_path="/Library/Fonts/ipaexg.ttf",
        width=600,
        height=300,
        prefer_horizontal=1,
        background_color='white',
        include_numbers=True,
        colormap='tab20',
        regexp=r"[\w']+",
    ).generate_from_text(tokenized_text)

# TF-IDFで Word Cloud作成
wc_1 = WordCloud(
        font_path="/Library/Fonts/ipaexg.ttf",
        width=600,
        height=300,
        prefer_horizontal=1,
        background_color='white',
        include_numbers=True,
        colormap='tab20',
    ).generate_from_frequencies(tfidf_dict)

# それぞれ可視化
plt.rcParams["font.family"] = "IPAexGothic"
fig = plt.figure(figsize=(12, 12), facecolor="w")
ax = fig.add_subplot(2, 1, 1, title="単語の出現回数で作成")
ax.imshow(wc_0)
ax.axis("off")
ax = fig.add_subplot(2, 1, 2, title="TF-IDFで作成")
ax.imshow(wc_1)
ax.axis("off")
plt.show()

出来上がった図がこちら。

明らかに下のやつの方がいいですね。