はじめての自然言語処理

前回はグラフベースのキーフレーズ抽出手法と pke での実験結果を紹介しました。今回は、spaCy, scattertext, ... 等々の OSS を用い各種のテキストマイニング手法についてコード例とサンプルプロットを交えながら説明したいと思います。

1. はじめに

本記事ではテキストマイニングの概要と代表的な手法について、コード例とサンプルプロットを交えて説明します。分析対象には、この連載で何度か用いている livedoor ニュースコーパスを用い、Google Colaboratory で動かすことを想定したコードスニペットを入れていきますので、実際に動かしたり対象を変えてみたりして試して頂けると良いかと思います。

2. テキストマイニングとは

テキストマイニングとは、ざっくり言うと「自然言語の文書データを対象に使用される単語の出現傾向等を分析して何らかの有益な情報を取り出すこと」となるでしょう。

「何らかの有益な」とは、ふわっとした表現ですが、テキストマイ二ングでは「何が有益なのか」を分析者が押さえていることが重要です。「何かを知りたい」や「何かを判断したい」というゴールがあり、その為の判断材料やヒントをテキストマイニングを用いて探し出します。ゴールに到達するには、どのような判断材料やヒントが必要か？、それらを得るには、どの分析手法が適切か？という思考で分析作業を進めていきましょう。

「とりあえずテキストマイニングツールを導入して、手元のデータを分析してみよう」とやると、いろいろな表やチャートが画面に表示され、「…それで？」となりがちです。漠然と「何かわかるんじゃないか？」ではなく、「どういうことを知りたいか／知れたらうれしいか」という意識があると分析結果の見え方も違ってくるのではないでしょうか。採掘（マイニング）するとき、自分が何を掘っているか分かっていないと折角の鉱物を見落としてしまうかもしれません。

適切な分析手法を選択したり分析結果を解釈したりするには、分析手法についてある程度の理解が必要です。数学的詳細までは必ずしも必要ないので、基本的な考え方や発想について次章以降で紹介していきます。

さて、少々前置きが長くなってしまったので、ここからは手を動かしなら進めていきましょう。

3. 分析環境のセットアップ

直前に、とりあえずツールを入れるのは如何なものか的なことを書いておいて何なのですが ¹ 、何はともあれ分析環境のセットアップをしていきましょう。とりあえず今回は spaCy（と GiNZA）をメインに可視化用ツール（word_cloud, scattertext）と日本語フォントを入れています。

以下、本記事中のコード例は Google Colaboratory でのコードセルの記述とします ² 。!で始まる行はコマンド実行、それ以外の行は Python3.x のコードです。

!pip install "https://github.com/megagonlabs/ginza/releases/download/v2.0.0/ginza-2.0.0.tar.gz"
!git clone https://github.com/amueller/word_cloud
!cd word_cloud && python setup.py install
!pip install scattertext jieba empath astropy gensim umap-learn
!apt install fonts-takao-pgothic

セットアップが終わったら一度ランタイムを再起動しておいてください。次は文書データを読みこみましょう。

4. データセットのロード

まずは分析対象の文書データとして livedoor ニュースコーパス ³ を取得して展開します。

!curl -O https://www.rondhuit.com/download/ldcc-20140209.tar.gz
!tar zxf ldcc-20140209.tar.gz

livedoor ニュースコーパスは9つのグループに分けられたニュース記事で構成されており、総記事数 7376 記事と比較的手ごろな大きさで扱いやすいコーパスです。 テキストマイニングでは文書データに付与された部署、日付、性別、商品カテゴリなどの属性データを用いて集計を行いますが、今回はニュース記事が属するグループを属性に見立てて分析をしていきます。

以下のようにして、展開した記事とグループをメモリ上にロードします。

import os
import numpy as np

LDCC_DIR = "./text"
groups = [f for f in os.listdir(LDCC_DIR) if os.path.isdir(LDCC_DIR + "/" + f)]

file_names = {}
for group in groups:
    files = os.listdir(LDCC_DIR + "/" + group)
    file_names[group] = files

texts  = []
labels = []

for group in file_names.keys():
    for file_name in file_names[group]:
        with open(os.path.join(LDCC_DIR, group, file_name), 'r') as f:
          lines = f.readlines()
          text = "".join(lines[3:])
        texts.append(text) 
        labels.append(group)
labels = np.array(labels) 

テキストマイニングではテキストを形態素解析で単語に分解し、その出現や係り受けをカウントします。 今回は spaCy/GiNZA を用いました。以下のようにしてtexts（文字列のリスト）を Docのリストにしてしまいましょう。Colaboratory だと 1 時間程度かかります。90分でセッションが切れるので直前で F5 を押してリロードしておく方が無難かもしれません。

import spacy
nlp = spacy.load('ja_ginza')
docs = list(nlp.pipe(texts, disable=['ner']))

5. 前準備

実際に分析を始める前にいくつか前準備をしておきます。

まず、GiNZA の形態素解析器である SudachiPy は “東京国際展示場” -> [“東京”, “国際”, “展示場”] のように比較的細かく名詞を分割してくるので、今回は以下のようにして名詞句を一つのトークンに再結合しました。これで“東京国際展示場"を 1 トークンとして扱えるようになります。

for doc in docs:
  with doc.retokenize() as retokenizer:
    for noun_chunk in doc.noun_chunks:
      retokenizer.merge(noun_chunk)  

次に単語の頻度をカウントして BOW ベクトル化するのに今回は sklearn の CountVectorizer を使うのでインポートしておきます。 その際、頻出語、希少語、ストップワードは除外しますが、これらの値は何度か使用するので変数にしておきます。

from sklearn.feature_extraction.text import CountVectorizer
NGRAM=1
MAX_DF=0.95
MIN_DF=0.03
NUM_VOCAB=1000

ストップワードは GiNZA が定義として持っているものに、本記事の執筆中でいくつか気になった単語を追加しています。

import ginza
stop_words = list(ginza.STOP_WORDS)
stop_words.extend(['max', 'エスマックス', 'smaxjp'])

今回は、「書か」、「書き」、「書け」等の活用形は「書く」の基本形(lemma)にして集計しています。 見た目の話ですが GiNZA の基本形は字面がやや固い印象⁴なので可視化の際に以下の関数を噛ましています。

def soften(word):
  replace_table = {
      '為る': 'する', '成る': 'なる', '遣る': 'やる', '有る': 'ある', '無い': 'ない',
      '御洒落': 'おしゃれ', '撫子': 'なでしこ', '未だ未だ': 'まだまだ', '迚も': 'とても',
      '唯': 'ただ', '筈': 'はず', '若し': 'もし'
  }
  return replace_table.get(word, word)

では、ようやくですが分析と可視化をしてきましょう。

6. 単語頻度

テキストマイニングの分析として一番単純なのは単語頻度でしょう。分析対象の中で沢山使われている単語は直感的にも重要な意味を持ちそうです。とりあえず、名詞を対象にして高頻度単語のトップ20を計算し、それらのグループ毎の内訳を見てみましょう。

CountVectorizer は NGRAM に対応していますが、今回は単語単体を扱うユニグラムで、ストップワードに指定した単語、95%以上の文書に出現する単語、3%以下の文書にしか出現しない単語は除外し、出現頻度の高い最大1000単語で BOW に変換します。

POS_NOUN = ['PROPN', 'NOUN'] # 固有名詞と名詞
tokens = []
for doc in docs:
    tokens.append(" ".join([token.lemma_ for token in doc if token.pos_ in POS_NOUN or len(POS_NOUN) ==0]))
cv = CountVectorizer(stop_words=stop_words, ngram_range=(1,NGRAM), max_df=MAX_DF, min_df=MIN_DF, max_features=NUM_VOCAB)
X_bow = cv.fit_transform(tokens).toarray()
print("Shape of X : %s" % (X_bow.shape,))
# Shape of X : (7376, 448)

結果として 7376文書が448次元の BOW に変換されました。BOW の次元と単語の対応は cv.vocabulary_ を使用します。 キーが単語、値が BOW のインデックスの辞書オブジェクトになっています。

vocab  = cv.vocabulary_ 
print("Num of vocab : %s" % (len(vocab)))
print("Sample of vocab : %s" % (list(vocab.keys())[:3]))
# Num of vocab : 448
# Sample of vocab : ['紹介', 'グーグル', 'トップ']

次は高頻度のトップ20単語のみのBOWベクトルに変換します。合わせてトップ20単語の文字列をリスト(words)にしておきましょう。

TOP_K = 20
sum_all = np.sum(X_bow, axis=0)
indices_topk = np.argsort(sum_all)[::-1][:TOP_K]
X_bow_topk = np.take(X_bow, indices_topk, axis=1)
print("Shape of X_bow_topk : %s" % (X_bow_topk.shape,))
# Shape of X_bow_topk : (7376, 20)

reverse_vocab = {vocab[k]:k for k in vocab.keys() }
words = [reverse_vocab[i] for i in indices_topk]

全体のトップ20単語のグループ毎の集計を計算します。

group_wordcount = {}
for group in groups:
    indices = np.where(labels == group)[0]
    X_group = np.take(X_bow_topk, indices, axis=0)
    sum_group = np.sum(X_group, axis=0) 
    group_wordcount[group] = sum_group

最後に以下のような関数でプロットです。

from matplotlib import pyplot as plt
import matplotlib.cm as cm
from matplotlib.font_manager import FontProperties
%matplotlib inline
plt.style.use('seaborn-whitegrid')

def plot_word_dist(words, group_wordcount, title, color="darkcyan"):
    left = np.array(list(range(len(words))))
    f, ax = plt.subplots(figsize=(6,3), dpi=180)
    fp = FontProperties(fname='/usr/share/fonts/truetype/fonts-japanese-gothic.ttf', size=6)
    bottom = np.zeros(len(words))
    cmap = cm.get_cmap('Set3')
    for i, group in enumerate(group_wordcount.keys()):
        vector = group_wordcount[group]
        color=1/len(groups) * i
        plt.bar(left, vector, label=group, bottom=bottom, color=cmap(color))
        bottom += vector
    ax.legend(loc='upper right')
    plt.xticks(left, words, rotation=290,  fontproperties=fp)
    plt.ylabel("出現頻度", fontproperties=fp)
    plt.title(title, fontproperties=fp)
    plt.tick_params(labelsize=6)
    plt.legend(fontsize=6)

plot_word_dist(words, group_wordcount, "コーパス全体の高頻度トップ20とニュースグループ毎の内訳")

最初は動詞、形容詞、形容動詞も含めてプロットしたのですが、トップ20単語のグループ内訳にあまり特徴がなかったので、名詞だけにしてみました。頻度を集計しただけですが、単語とグループの関係性をさらっと俯瞰できますね。

次は、ある部分集合に特徴的な単語を抽出してみましょう。

7. JLHスコア

JLHスコアはコーパスの部分集合における重要度を示す指標です。時系列になっているコーパスであれば「最近1週間での話題のワード」のような使い方ができるでしょう。計算は比較的簡単で以下のとおりです。

ここで、絶対割合変化と相対割合変化は以下のとおりです。

• 絶対割合変化 = 部分集合の出現率 - コーパス全体の出現率
• 相対割合変化 = 部分集合の出現率 / コーパス全体の出現率

今度は名詞、動詞、形容詞、形容動詞で計算してみましょう。

POS = ['PROPN', 'NOUN', 'VERB', 'ADJ', 'ADV']
tokens = []
for doc in docs:
    tokens.append(" ".join([token.lemma_ for token in doc if token.pos_ in POS or len(POS) ==0]))
cv = CountVectorizer(stop_words=stop_words, ngram_range=(1,NGRAM), max_df=MAX_DF, min_df=MIN_DF, max_features=NUM_VOCAB)
X_bow = cv.fit_transform(tokens).toarray()
print("Shape of X : %s" % (X_bow.shape,))
#Shape of X : (7376, 766)
vocab  = cv.vocabulary_ 

コーパス全体での単語の出現率を計算します。

X_sum = np.sum(X_bow, axis=0)
total_occurence_rate = X_sum / len(X_bow)

ニュースグループ「独女通信」の単語の出現率を計算します。

group = "dokujo-tsushin"
indices = np.where(labels == group)[0]
count = len(indices)
X_group = np.take(X_bow, indices, axis=0)
group_sum = np.sum(X_group, axis=0)
group_occurence_rate = group_sum / count

JLH スコアは先ほどの数式の通り、以下のようにして計算します。

sub = group_occurence_rate - total_occurence_rate
div = group_occurence_rate / total_occurence_rate
jlh = sub * div

棒グラフばかりでは面白くないので（？）、今度は JLH スコアのトップ20をワードクラウドにして見ましょう。 描画ロジックはこんな感じです。

from wordcloud import WordCloud
def plot_wordcloud(vocab, jlh_indices_topk, jlh_score_topk, group):
  reverse_vocab = {vocab[k]:k for k in vocab.keys() }
  words_topk = [soften(reverse_vocab[i]) for i in jlh_indices_topk]

  word_freq = {}
  for i, word in enumerate(words_topk):
    word_freq[word] = jlh_score_topk[i]

  fpath = '/usr/share/fonts/truetype/fonts-japanese-gothic.ttf'
  fp = FontProperties(fname=fpath, size=16)
  wordcloud = WordCloud(background_color="white",font_path=fpath, width=900, height=500)
  image = wordcloud.generate_from_frequencies(word_freq)
  plt.figure(figsize=(15,12))
  plt.imshow(image)
  plt.axis("off")
  plt.title("グループ %s における JLH スコア トップ20" % group, fontproperties=fp)
  plt.show()

TOP_K_JLH = 20
jlh_indices_topk = np.argsort(jlh)[::-1][:TOP_K_JLH]
jlh_topk = np.take(jlh, jlh_indices_topk)
plot_wordcloud(vocab, jlh_indices_topk, jlh_topk, group)

確かにそれっぽい結果ですね。計算間違いはしてなさそうです。

正直、ワードクラウドは見た目が派手なだけで分かりにくいので、あまり好きではありません。 JLH のようにコーパスの部分集合に特徴的な単語を示す指標として Scaled F-Score とその描画 ツールがありますので、そちらの紹介に移りましょう。

8. Scaled F-Score

Scaled F-Score はカテゴリと高い関連性のある単語をカテゴリ固有適合率とカテゴリ固有生起頻度の調和平均で評価する手法です。 単語 i 、カテゴリ j のカテゴリ固有適合率を prec(i, j) 、カテゴリ固有生起頻度を freq(i, j) 、単語のカテゴリ内での出現回数を #(w_i, c_j) とすると、Scaled F-Score H_β(i, j) は以下のようになります。

ようするに着目する単語、カテゴリがあり、以下の二つの要素のバランスを取って評価しているということです。

• 単語とカテゴリとの関連度：着目単語が着目カテゴリ内で他カテゴリに比べどの程度出現しているか
• カテゴリ内での単語の重要度：着目単語が着目カテゴリ内で他単語に比べてどの程度出現しているか

それでは実際にプロットしてみましょう。 scattertext ⁵を用いるとラベル文字列が重ならないように調整された美しくインタラクティブな散布図を簡単に描画できます。scattertext は spaCy を利用していますが、GiNZA があれば大丈夫です。

まずはニュースグループ「スポーツウォッチ」と「独女通信」のデータを pandas の DataFrame にします。

import pandas as pd
import scattertext as st
ldcc_df = pd.DataFrame([[example[0], example[1]] for example in zip(labels,docs) 
    if example[0] in ["dokujo-tsushin", "sports-watch"]])
ldcc_df.columns = ["group", "text"]

DataFrameの中身はこんな感じです。今回は全文書がDocに変換済みだったのでtext列をDoc型にしています ⁶。

print(ldcc_df.head(3))
#            group                                               text
#0  dokujo-tsushin  (新しい, コミュニケーション, の, ひと, つ, の, 形, と, し, て, 、, 今...
#1  dokujo-tsushin  (昨年, 12, 月, の, 時点, で, 、, 日本, に, おけ, る, SNS, の,...
#2  dokujo-tsushin  (人, は, 色々, な, 瞬間, に, 絶望, を, 感じる, 。, 仕事, で, 失敗,...

次は DataFrame を st.CorpusFromParsedDocuments() で ParsedCorpus に変換します。このとき spaCy の Doc から特徴量を抽出するクラスを feats_from_spacy_doc パラメータで指定できます。 デフォルトの scattertext.features.FeatsFromSpacyDoc を使用しても良いのですが、品詞でフィルタを掛けたかったので以下のようなクラスを作りました。

from collections import Counter
from itertools import chain

class PosFilteredUnigramFeatsFromSpacyDoc(object):
    def __init__(self,
                 use_lemmas=False,
                 strip_final_period=False,
               poses_to_include = ['PROPN', 'NOUN', 'VERB', 'ADJ', 'ADV']):
        self._use_lemmas = use_lemmas
        self._strip_final_period = strip_final_period
        self._poses_to_include = poses_to_include

    def _post_process_term(self, term):
        if self._strip_final_period and (term.strip().endswith('.') or term.strip().endswith(',')):
            term = term.strip()[:-1]
        return term

    def get_doc_metadata(self, doc):
        return Counter()

    def get_feats(self, doc):
        ngram_counter = Counter()
        for sent in doc.sents:
            unigrams = self._get_unigram_feats(sent)
            ngram_counter += Counter(chain(unigrams, []))
        return ngram_counter

    def _get_unigram_feats(self, sent):
        unigrams = []
        for tok in sent:
            if tok.pos_ in self._poses_to_include:
                if self._use_lemmas and tok.lemma_.strip():
                    unigrams.append(self._post_process_term(soften(tok.lemma_.strip())))
                elif tok.lower_.strip():
                    unigrams.append(self._post_process_term(tok.lower_.strip()))
        return unigrams

    def has_metadata_term_list(self):
        return False

    def get_top_model_term_lists(self):
        raise Exception("No topic models associated with these features.")

後は以下のようにしてParsedCorpus に変換し,

import scattertext as st
from scattertext.features import FeatsFromSpacyDoc
feats_from_spacy_doc=PosFilteredUnigramFeatsFromSpacyDoc(use_lemmas=True)
#feats_from_spacy_doc=FeatsFromSpacyDoc() # こちらは scattertext のデフォルト実装
corpus = st.CorpusFromParsedDocuments(ldcc_df, category_col="group", parsed_col="text", 
                                      feats_from_spacy_doc=feats_from_spacy_doc).build()

HTML にして、

html = st.produce_scattertext_explorer(corpus, category='sports-watch',
          category_name='スポーツウォッチ', not_category_name='独女通信',
          width_in_pixels=1000, minimum_term_frequency=50, max_terms=4000)

描画します⁷。散布図の座標はそれぞれのカテゴリにおける単語の頻度、青～黄～赤の色が Scaled F-Score に対応します。ドットの密度が濃い部分はラベルが表示されていませんがマウスカーソルを合わせるとポップアップ表示されます。

from IPython.display import HTML
HTML(html)

ここからは TV のニュースなどでよく見かける単語共起ネットワークを描画してみましょう。

9. 単語共起ネットワーク

単語共起ネットワークは文や段落において同時に出現（共起）する単語ペアをネットワーク状に表示して文書集合の傾向を可視化するものです。TV のニュースで SNS 上のつぶやきを分析したものを目にしたことがあるかと思います。

どこまでの範囲で同時に出現したら共起とするかは、文や段落とすることが多いようです。今回は文単位で計算します。 まずは、文単位に” “で分かち書きしたリスト（とラベル）に変換します。

tokens = []
labels_sent = []
for i, doc in enumerate(docs):
    for sent in doc.sents:
        tokens_of_sent = [token.lemma_ for token in sent if token.pos_ in POS or len(POS) ==0]
        if len(tokens_of_sent) < 2:
            continue
        tokens.append(" ".join(tokens_of_sent))
        labels_sent.append(labels[i])
labels_sent = np.array(labels_sent)

次に文中での単語共起ペアの計算ですが、今回は文単位の BOW の outer で以下のように計算します。BOW の長さが n なら、n x n の2次元配列で単語 i と単語 j が共起していれば回数に関係なく (i, j) = 1、共起がなければ 0 としました。

def calc_cooccurrence(bow):
    length = bow.shape[0]
    clear_diag = np.ones(length, dtype=np.int) - np.eye(length, dtype=np.int)
    array = np.minimum.outer(bow, bow) * clear_diag
    return (array > 0).astype(int)

それではニュースグループ「ムービーエンタ」の単語共起ネットワークを表示してみましょう。ポイントは二つあります。

一つは共起の指標に共起回数ではなく、Jaccard係数を使うことです。 Jaccard係数は前回出てきたので詳しい説明は省きますが、単語A, B の共起であれば、単語AもしくはBを含む文の集合のうち、AとBの両方を含む文の割合になります。

例えば、文数100、単語 A、B、C があるとして (A, B) 、(A, C)の共起が共に10個の文で出現したとします。このとき単語 A、B、C はそれぞれ10, 95, 11個の文に出現していたとするとどうでしょう。共起の回数自体は (A, B) も (A, C) も10回づつで同じですが、「対で出現する傾向の強さ」は同じとは言い難いことが直感的にわかると思います。そこで Jaccard係数で回数ではなく割合に着目する訳です。

もう一つのポイントは低頻度の単語を強めにフィルタしておくことです。 Jaccard係数で頻度ではなく割合を見ることにしました。先の例の続きで文数100の中で低頻度語 X, Y が1回づつ、同一文の中に出現するケースでは、Jaccard係数 = 1.0 の最大値になってしまいます。これでは描画内容が低頻度語で散らかってしまうので、予めフィルタで落としておくわけです。 以下では min_df で文（記事ではなくて）に出現する確率が1%以下の単語を切り捨てています。

group = "movie-enter"
indices = np.where(labels_sent == group)[0]
tokens_group = np.take(tokens, indices, axis=0)
cv_group = CountVectorizer(stop_words=stop_words, ngram_range=(1,NGRAM), min_df=len(indices) // 100, max_features=NUM_VOCAB) 
X_bow_group = cv_group.fit_transform(tokens_group).toarray()
vocab_group  = cv_group.vocabulary_ 
vocab_size = X_bow_group.shape[1]
intersection = np.zeros([vocab_size, vocab_size]).astype(int)
for bow in X_bow_group:
    intersection += calc_cooccurrence(bow)
sum_occurrences = np.sum((X_bow_group > 0).astype(int), axis=0)
union = np.add.outer(sum_occurrences, sum_occurrences) - intersection
union[np.where(union==0)[0]] = 1 # to avoid division by zero.
jaccards = intersection / union

Jaccard係数が計算できたので、次は NetworkX を使ってグラフを構築する関数です。

import networkx as nx
from itertools import combinations
def create_graph(scores, vocab, topk=50):
    vocab = {soften(word):id for word,id in vocab.items()}
    G = nx.Graph()
    G.add_nodes_from([w for w in vocab.keys()])

    flatten_scores = scores.reshape([-1])
    # "*2" は scores は同じ値が対角線両側の二回ずつ入る為
    thresh = np.sort(flatten_scores)[::-1][topk*2] 

    reverse_vocab = {vocab[k]:k for k in vocab.keys() }

    # Top K 以上の score(Jaccard係数) の単語ペアにエッジを張る
    for (i,j) in combinations(range(scores.shape[0]), 2):
        weight = scores[i][j]
        if weight > thresh:
            G.add_edge(reverse_vocab[i], reverse_vocab[j], weight=weight)

    # エッジが張られていない頂点を削除
    isolated = [n for n in G.nodes if len([i for i in nx.all_neighbors(G, n)]) == 0]
    for n in isolated:
       G.remove_node(n)
    return G

構築したグラフを使ってネットワーク図を描く関数です。ネットワーク中のエッジの太さは Jaccard係数を、 ノードの大きさには Jaccard係数をエッジの重みとして計算した PageRank の値を反映させています。 エッジで接続された一固まり毎に色を変えると見栄えが良くなります。

import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from matplotlib.font_manager import FontProperties
%matplotlib inline

def plot_network(group, G, k=None, iterations=50, fontsize=12, node_size_factor=30000, edge_width_factor=20):
    fp = FontProperties(fname='/usr/share/fonts/truetype/fonts-japanese-gothic.ttf', size=fontsize)
    plt.figure(figsize=(12, 8))

    layout = nx.spring_layout(G, k=k, iterations=iterations)  

    pr = nx.pagerank(G)
    pr_values = np.array([pr[node] for node in G.nodes()])

    connecteds = []
    colors = []
    for i, c in enumerate(nx.connected_components(G)):
      connecteds.append(c)
      colors.append(1/18 * i)

    node_colors = []
    for node in G.nodes():
      for i, c in enumerate(connecteds):
        if node in c:
          node_colors.append(colors[i])
          break

    nx.draw_networkx_nodes(G, layout, node_color=node_colors, cmap=plt.cm.get_cmap('Set3'), 
        alpha=0.7, node_size=pr_values * node_size_factor)

    labels = nx.draw_networkx_labels(G, layout)
    for t in labels.values():
        t.set_fontproperties(fp)

    edge_width = [d["weight"] * edge_width_factor for (u, v, d) in G.edges(data=True)]
    nx.draw_networkx_edges(G, layout, alpha=0.4, edge_color="darkgrey", width=edge_width)
    plt.title("グループ %s における単語共起ネットワーク" % group, fontproperties=fp)
    plt.axis('off')

ようやくですが、こんな感じの表示になります。k や iterations は NetworkX の spring_layout のパラメータです。 描画結果を見ながら適当に決めてます。

import math
G = create_graph(jaccards, vocab_group, topk=30)
k = 2.0/math.sqrt(len(G.nodes()))
plot_network(group, G, k=k, iterations=65)

ここまでは文書に予め付与された属性（今回の場合はニュース記事のグループ）で対象を絞り込んだり、集計したりしましたが、ここからは分析によって新たな属性を見出すということをして見ましょう。

10. クラスタリング

クラスタリングは分析対象の集合を類似度などに基づいて複数の部分集合（＝クラスタ）に分割する手法で、未知のコーパスの全体概要を把握したい場合などに有用です。

分割はパラメータを指定して後はアルゴリズム任せですので、分割された結果はクラスタ1、クラスタ2、… でしかなく「スポーツ」、「経済」のような名前があるわけではありません。前述の JLH スコア等で各クラスタの特徴語を確認して人手で名前を付けるなりなんなりすることになります。後はその分類先のクラスタを属性として、内訳を見るなり散布図を作るなり、ここまで述べてきた手法で分析しましょう。

クラスタリングの手法には大別して階層型と非階層型に分けられ、それぞれ以下のような特徴があります。

• 階層型(群平均法、ウォード法)
  • 距離の近いクラスタ併合し、階層的に併合を繰り返す。
  • 併合過程を確認できるが、大量データの場合に計算量が問題になる。
• 非階層型（K-Means法）：
  • 大量データに強いが、最初に分類数を指定する必要がある。

階層型については前回少し触れたので、今回は非階層型手法の K-Means法を利用してみましょう。

K-Means法は最初に分割数(k)を指定します。後は以下のように最初は k 個の重心をランダムに決定（①）、全ての点を最も近い重心に紐付け k 個のクラスタを形成（②）、クラスタの形成結果から新に k 個の重心を決定（③）、全ての点を最も近い重心に紐付け（④）…という具合に重心が移動しなくなるまで続けます。

手法の概要が分かったところで、実際にクラスタリングをしてみましょう。重心や距離を計算するためにはコーパスの文章をなんらかのベクトルにしなければなりませんが、ここは本連載の第1回で紹介した TF-IDF でベクトル化してみます。

from sklearn.feature_extraction.text import TfidfTransformer
tokens = []
for doc in docs:
    tokens.append(" ".join([token.lemma_ for token in doc if token.pos_ in POS or len(POS) ==0]))
cv = CountVectorizer(stop_words=stop_words, ngram_range=(1,NGRAM), max_df=MAX_DF, min_df=MIN_DF, max_features=NUM_VOCAB)
X_bow = cv.fit_transform(tokens).toarray()
vocab  = cv.vocabulary_ 
X_tfidf = TfidfTransformer().fit_transform(X_bow).toarray()
print("Shape of X_tfidf : %s" % (X_tfidf.shape,))
# Shape of X_tfidf : (7376, 766)

文章を TF-IDF ベクトルにしたら K-Means法でクラスタリングします。今回は sklearn の実装を使いました。 分割数はコーパスのニュースグループの数に合わせて 9、残りのパラメータはとりあえずデフォルトのままです。

from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=9, random_state=0).fit(X_tfidf)

クラスタリングの結果は以下のようにして取り出します。

cluster_labels = kmeans.labels_
print("Shape of cluster_labels : %s" % (cluster_labels.shape,))
# Shape of cluster_labels : (7376,)

ただし、K-Means法は最初の初期値の配置によって結果が変わってきます。何度か実行しクラスタ内誤差平方和（SSE :Sum of Squared Errors）の小さいもの（＝クラスタが良くまとまっているもの）を選ぶと良いでしょう。クラスタ内誤差平方和はクラスタ内の全点と重心の平方距離の和で、sklearn を使えば以下のようにして確認できます。

print("SSE = %.2f" % (kmeans.inertia_))
# SSE = 6136.37

せっかく、クラスタリングしたので PCA で次元削減してプロットしてみましょう。

from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties

def draw_scatter_plot(vectors, cluster_labels, label_names=None, 
                      with_index=False, xlim=None, ylim=None, title = "", aspect=(5,4), dpi=150, marker_size=2):
    fp = FontProperties(fname='/usr/share/fonts/truetype/fonts-japanese-gothic.ttf', size=8)
    pca = PCA(n_components=2)
    coords = pca.fit_transform(vectors)
    print("shape of coords : %s" % (coords.shape,))
    fig, ax = plt.subplots(figsize=aspect, dpi=dpi)
    for c in np.unique(cluster_labels):
        indices = np.where(cluster_labels == c)[0]
        coords_c = np.take(coords, indices, axis=0)
        if label_names :
          label = label_names[c]
        else:
          label = "cluster:%d" % c     
        ax.scatter(coords_c[:, 0], coords_c[:, 1], s=marker_size, label=label) 
    ax.legend(loc='upper right', fontsize=6)
    plt.title(title, fontproperties=fp)
    if xlim:
        plt.xlim(xlim)
    if ylim:
        plt.ylim(ylim)
    plt.tick_params(axis='y', which='major', labelsize=5)    
    plt.tick_params(axis='x', which='major', labelsize=5)
    plt.show()

 draw_scatter_plot(X_tfidf, cluster_labels)  # クラスタリング結果のプロット

 group_labels = [groups.index(label) for label in labels]
 draw_scatter_plot(X_tfidf, group_labels, groups)  # 元々のグループ分けのプロット

左側がクラスタリング結果のプロットで、右側が元々のグループ分けに基づいて色分けした結果です。左と右で色の関係はありません。左のピンクと右の赤、左の紫と右のグレーを見ると、そこそこ元のグループ分けに似通ったクラスタが形成されているようです。

ところで、今回は元々のニュースグループの数が 9 だったので、さらっと分類数を 9 にしましたが、本当に未知のコーパスを分析するときはどうやって分割数を決定したらよいでしょう？ 一般的に良く紹介されているのはエルボー法という手法です。話はシンプルで分割数を変えながら複数回 K-Means法を実行し、分割数と SSE の関係をプロットします。ブロットの折れ線を腕に見立て肘（エルボー）のように折れ曲がっている部分を選ぶのが良いとされています。

以下は今回のデータでエルボー法を試したプロットです。分割数 9 のところが肘といえば肘でしょうか。ちなみにユリウス・カエサルは『人間はみな自分の見たいものしか見ようとしない。』と言ったそうです。。。現実のデータで試すと分かりやすい肘が現れることは、なかなかないらしいです。

クラスタリングではコーパスを似た者同士の部分集合に分割する処理でしたが、次は文章を構成する話題を推定するトピック分析を試してみましょう。

11. トピック分析

トピック分析は文章がどのような話題（トピック）で構成されているかを推定するものです。LDA (Latent Dirichlet Allocation) はその一手法で、文書・トピック・単語の関係を確率的に扱う言語モデルになります。理論的なところはさらっとは説明できないので省きますが、 LDA による文書生成は以下のようなイメージになります。

まず、文書毎にガラポンを回してトピックを決める為のサイコロを得ます。ガラポンから出てくるサイコロは毎回形が異なり、振った時に出る目（トピック）の確率も違います（①）。次に①のサイコロを振って出た目（トピック）に対応したサイコロを選びます。このサイコロはトピック毎に用意されており、それぞれ出る目（単語）の確率が異なります（②）。そして②のサイコロを振って単語を選びます（③）。文書の単語数だけ②、③を繰り返します。文書が変われば再度ガラポンでトピックを決めるサイコロを選びなおします（④）。

LDA での学習とは「ガラポンを回したときにどのような形のサイコロが出てくるか」、「トピック毎に用意されたサイコロはそれぞれどのような形をしているか」という具合を、文書の生成結果が学習対象のコーパスに出来るだけ近づくよう調整するということになります。

学習コーパスに対し、トピック数を指定して学習し LDA のモデルが出来上がると、以下のようなことがわかるようになります。

• コーパスから潜在的なトピックが抽出できる。
• 各トピックにおける単語の出現確率がわかる。
• 文書に対してトピックの混合比が推論できる。

さて、トピック数としていくつを指定するかですが、 LDA の性能指標としては以下の二つがあります。

Perplexity

Perplexity は「モデルの困惑度合い」を示す指標です。言語モデルが文書を生成するとき語彙中から単語を選ぶことになります。この単語選択の確信度の逆数が Perplexity になります。確信度が 0.5 (=½) であれば、 Perplexity = 2 という訳です。数式で示すと以下のようになります⁸。

• D_test : テストデータの文書集合
• M : D_test の文書数
• w_d : 文書 d の単語集合
• N_d : 文書 d の単語数

Perplexity は言語モデルとしての予測性能を見ているわけですが、Perplexityに基づいてトピック数を選ぶと人の感覚とは異なる結果になるという報告もあるようです⁹。 直接的に推定されたトピックの質を評価する手法として Coherence がありますので、そちらも見てみましょう。

Coherence

Coherence は簡単に言ってしまえば、トピックを代表する単語の一貫性度合いになります。例えば、{"ザキトワ"、"ダービー”, “全英オープン”} であれば、何やらスポーツに関するトピックのようです。しかし、{“博多駅"、"年末年始"、"ビッグバード”} では何のことだか良くわかりません。前者は Coherence が高く、後者は低いということになります。

Coherence には様々な種類があるのですが、今回は学習コーパス以外のデータを用意する必要がなく、高速、手軽な U Mass Coherence を使用してみます¹⁰。

• M : 各トピックを代表する単語数
• V^{(t) = ( v₁(t), v₂(t), … v_M(t)) : トピック t を代表する M 個の単語}
• D(v₁, v₂): v₁, v₂ が共起する文書数　

それではこちらも実際に動かしてみましょう。

livedoor ニュースコーパスでの実験

まずは LDA の入力データとして TF-IDF ベクトルを作ります。本来 LDA は BOW で入力するもののようなのですが、結果を目視で確認したところ、 TF-IDF を用いたほうがやや良質のトピックを抽出できていたので、今回はそのようにしています。

from sklearn.feature_extraction.text import TfidfVectorizer
tokens = []
for doc in docs:
    tokens.append(" ".join([token.lemma_ for token in doc if token.pos_ in POS or len(POS) ==0]))
tv = TfidfVectorizer(stop_words=stop_words, ngram_range=(1,NGRAM), max_df=MAX_DF, min_df=MIN_DF, max_features=NUM_VOCAB)
X_tfidf = tv.fit_transform(tokens)
vocab = tv.vocabulary_
reverse_vocab = { vocab[k]:k for k in vocab.keys() }
vocab_size = len(vocab)

トピック数の指定ですが、とりあえず 5 を指定しました。max_iter はデフォルトで 10 でしたが少し足らないようでしたので 30 を指定して、 Perplexity の下落幅はデフォルトのままで 0.1 を切ったら停止です¹¹。

from sklearn.decomposition import LatentDirichletAllocation as LDA
lda = LDA(n_components=5, max_iter=30, n_jobs=-1, verbose=1, evaluate_every=1)
lda.fit(X_tfidf) 

# iteration: 1 of max_iter: 30, perplexity: 1098.6809
# iteration: 2 of max_iter: 30, perplexity: 931.0230
# ...
# iteration: 20 of max_iter: 30, perplexity: 841.1229
# iteration: 21 of max_iter: 30, perplexity: 840.9431
# iteration: 22 of max_iter: 30, perplexity: 840.8543
# LatentDirichletAllocation(batch_size=128, doc_topic_prior=None,
#                          evaluate_every=1, learning_decay=0.7,
#                          learning_method='batch', learning_offset=10.0,
#                          max_doc_update_iter=100, max_iter=30,
#                          mean_change_tol=0.001, n_components=5, n_jobs=-1,
#                          perp_tol=0.1, random_state=None,
#                          topic_word_prior=None, total_samples=1000000.0,
#                          verbose=1)

U Mass Coherence の値も見てみましょう。

U Mass Coherence の計算はこんな感じです。ちなみに今回は LDA の実装に sklearn を使ったので、自前で実装しましたが、 gensim には U Mass Coherence を含む複数の実装が組み込まれています。

def calc_cooccurrence(bow):
    length = bow.shape[0]
    clear_diag = np.ones(length, dtype=np.int) - np.eye(length, dtype=np.int)
    array = np.minimum.outer(bow, bow) * clear_diag
    return (array > 0).astype(int)

import math
def calc_umass_coherence(v_t, d_vl, d_vm_vl, as_mean=False):
    sum_of_word_pairs = 0.0
    num_of_word_pair = 0
    for m in range(1, len(v_t)):
        for l in range(m):
            sum_of_word_pairs += math.log((d_vm_vl[v_t[m]][v_t[l]] + 1.0) / d_vl[v_t[l]])
            num_of_word_pair +=1
    if as_mean:
      return sum_of_word_pairs / num_of_word_pair
    else:
      return sum_of_word_pairs

推定したトピック毎にトップ10単語を抽出し、各トピックの U Mass Coherence を計算して平均しました。

n_components = 5
X = X_tfidf.toarray()
d_vm_vl = np.zeros([vocab_size, vocab_size]).astype(int)
for bow in X:
  d_vm_vl += calc_cooccurrence(bow)
d_vl = np.sum((X > 0).astype(int), axis=0)
sum_coherences = 0.0
for i, component in enumerate(lda.components_):
  top_k_indices = component.argsort()[::-1][:10]
  sum_coherences += calc_umass_coherence(top_k_indices, d_vl=d_vl, d_vm_vl=d_vm_vl)
u_mass = sum_coherences / n_components
u_mass
# -46.755651246363875

数字だけ見てもよくわからないのでプロットしてみましょう。まずはトピック毎にトップ10単語を抽出します。

top_k_words = []
for i, component in enumerate(lda.components_):
    top_k_indices = component.argsort()[::-1][:10]
    top_k_word = {soften(reverse_vocab[k]): component[k] for k in top_k_indices}
    top_k_words.append(top_k_word)

各トピックは全語彙に対する確率分布を持っていますので、そのベクトルを T-SNE で次元削減して散布図としてプロットします。ただ色付きの点をプロットしてもつまらないので、先ほどのトップ10単語を円状のワードクラウドにして配置してみました。

from matplotlib import pyplot as plt
import matplotlib.cm as cm
from matplotlib.font_manager import FontProperties
%matplotlib inline
plt.style.use('seaborn-whitegrid')
from sklearn.manifold import TSNE
from wordcloud import WordCloud
import random
from functools import partial

def topic_cloud(word_freqs, hue_index, size=300):  

  def circle(size=300):
    d = size
    r = size // 2
    m = r * 0.86
    x, y = np.ogrid[:d, :d]
    mask = (x - r) ** 2 + (y - r) ** 2 > m ** 2
    return 255 * mask.astype(int)

  def random_color_func(word, font_size, position, orientation, 
                        random_state=None, num_topics=9, hue_index=0, **kwargs):
    fluc = 30
    hue = 360 // num_topics * (hue_index) + fluc
    return "hsl(%d, 100%%, %d%%)" % (hue + random.randint(-fluc, fluc), random.randint(30, 40))

  fpath = '/usr/share/fonts/truetype/fonts-japanese-gothic.ttf'
  fp = FontProperties(fname=fpath, size=16)
  wordcloud = WordCloud(background_color="white",font_path=fpath, mask=circle(size))
  wc = wordcloud.generate_from_frequencies(word_freqs[hue_index])
  color_func = partial(random_color_func, num_topics=len(word_freqs)  ,hue_index=hue_index)
  return wc.recolor(color_func=color_func)

def draw_lda_cloud(components, top_k_words, aspect=(10,8), dpi=180, cloud_size_pct=0.2):

  def draw_topic_cloud(x, y, top_k_words, index, size=100):
    tc = topic_cloud(top_k_words, hue_index=index, size=300)
    r = size//2
    plt.imshow(tc, extent=(x - r, x + r, y -r, y + r))

  tsne = TSNE(n_components=2)
  X_2d = tsne.fit_transform(components)

  xmin = X_2d[:,0].min()
  xmax = X_2d[:,0].max()
  ymin = X_2d[:,1].min()
  ymax = X_2d[:,1].max()
  xmargine  = (xmax - xmin) * 0.3
  ymargine = (ymax - ymin) * 0.3
  cloud_size = ((xmax - xmin) + (ymax - ymin)) / 2 * cloud_size_pct

  fig, ax = plt.subplots(figsize=aspect, dpi=dpi)

  for i in range(len(X_2d)):
    #ax.scatter(X_2d[i][0], X_2d[i][1], s=4, label=i) 
    draw_topic_cloud(X_2d[i][0], X_2d[i][1], top_k_words, i, size=cloud_size)

  plt.xlim((int(xmin - xmargine), int(xmax + xmargine)))
  plt.ylim((int(ymin - ymargine), int(ymax + ymargine)))
  #plt.axis("off")
  plt.tick_params(grid_alpha=0.6, grid_linewidth=0.5)
  plt.show()

 draw_lda_cloud(lda.components_, top_k_words, cloud_size_pct=0.4)

ざっと見ると、独女通信、スポーツウォッチ、ムービーエンタっぽいトピックが形成され、残りの二つのトピックはその他のグループの混ぜ物になっているようです。このプロットでは見た目で分かる以上のことはわかりませんが、可視化ツールの pyLDAvis を用いるとさらに詳しく確認することができます。 pip でインストールし、

!pip install pyLDAvis

以下のようにして表示します。日本語でも特に問題は起こりませんでした。

import pyLDAvis
import pyLDAvis.sklearn
pyLDAvis.enable_notebook()
panel = pyLDAvis.sklearn.prepare(lda, X_tfidf, tv, mds='tsne')
panel

以下はスクリーンショットですが、実際はトピックを表す左側の円にマウスを合わせると、選択したトピックに紐つく単語が右側に表示されますので、LDA でトピック分析を行う場合は合わせて利用してみてはいかがでしょうか。

先ほどは、トピック数をいきなり 5 と決めましたが、実はコーパスを 4:1 で学習データとテストデータに分け、トピック数を変更しながら、Perplexity と U Mass Coherence を計測してみました。各トピック数で5回試行し、その平均値を使っています。Perplexity は小さい方が良い指標、Coherence は大きい方が良い指標になります。

上記は LDA への入力に TF-IDF ベクトルを入れた結果です。Perplexity （というか TF-IDF で Perplexity 的な計算をした値）は右肩上がりですが、入力に BOW を用いると右肩下がりの傾向になりました。

じつは U Mass Coherence 的には 3 がベストだったのですが、トピック数 3 だと分析的につまらない気がしたので 5 を選んでしまいました。テキストマイニングをするときは正解が明確にあるわけでもないので、目安として使えば良いのではないかと思います。

12. 類義語の自動抽出

最後に類義語の自動抽出もやってみましょう。と言ってもコーパスから単語ベクトルを学習し、それに基づいて階層的クラスタリングで組み上げる過程を可視化するだけですので、そんな大した話でもありません。

まずは Doc を分かち書きテキストとして出力し、

with open("ldcc_wakati.txt", "w") as f:
  for doc in docs:
    tokens = [token.text for token in doc]
    tokens.append("\n")
    f.write(" ".join(tokens))

fastText で単語ベクトルの学習を行います。

!pip install gensim
!git clone https://github.com/facebookresearch/fastText
!cd fastText && make
!./fastText/fasttext skipgram -input ./ldcc_wakati.txt -output fasttext_ldcc -epoch 10 -dim 100 -minCount 20
# Read 4M words
# Number of words:  11020
# Number of labels: 0
# Progress: 100.0% words/sec/thread:   36471 lr:  0.000000 avg.loss:  2.028130 ETA:   0h 0m 0s

単語ベクトルが出来上がったら、gensim でロードして、

from gensim.models import KeyedVectors
wv = KeyedVectors.load_word2vec_format('./fasttext_ldcc.vec', binary=False) 

以下のようにして描画します。

import pandas as pd
from matplotlib import pyplot as plt
import matplotlib.cm as cm
from matplotlib.font_manager import FontProperties
%matplotlib inline
plt.style.use('seaborn-whitegrid')

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage, dendrogram

def draw_similar_word_dendrogram(wv, word, topn=30, method='average',  metric='cosine'):
    words = [x[0] for x in sorted(wv.most_similar(word, topn=topn))]
    words.append(word)
    vectors = [wv[word] for word in words]
    df = pd.DataFrame(vectors, index=words)
    lm = linkage(pdist(df, metric=metric), method=method)

    f, ax = plt.subplots(figsize=(5,3), dpi=140)
    fp = FontProperties(fname='/usr/share/fonts/truetype/fonts-japanese-gothic.ttf', size=8)
    dendrogram(lm, labels=words, orientation='right', ax=ax, color_threshold=0.0)
    for label in plt.gca().get_yticklabels():
        label.set_fontproperties(fp)

draw_similar_word_dendrogram(wv, "スマホ", 10)

「スマホ」と類似度が高いトップ10単語の樹形図を表示すると以下のようになりました。

13. おわりに

今回は少し趣向を変え、テキストマイニングで用いられる手法を OSS を用いてあれこれ試してみました。アドホックに分析するときは GUI のあるツールの方が便利ですし、フリーで使えるツールもあります。ただ繰り返してパラメータを探索する場合などプログラムで記述出来る方が便利な場合もありますので、そこは使い分けかと思います。次回は感情極性辞書を用いた感情分析についてご紹介する予定です。