1. 本Part概要

前PartではB'zの歌詞を「TF-IDF」を用いた分析を行いました。
本Partではトピックモデルの一つである「LDA」を用いた分析についてお話しします。

pira-nino.hatenablog.com

2. LDAとは

2.1 LDAのイメージ

先に簡単な説明をしてしまいます。 LDAは「たくさんの文書データから単語のグルーピングを行う」モデルです。
このグループ1つ1つを「トピック」と呼びます。

例えば、大量のニュース記事にLDAを適用する例を考えます。

f:id:pira_nino:20180729204436p:plain — ニュース記事データにLDAを適用した例

LDAでは「各トピック（トピック数は予め指定）における各単語の所属確率」が算出されます。
理論的なことはさておき、文書データから単語をいくつかのグループに自動で分けてくれる手法 との理解で大丈夫です。

よく勘違いされることとして以下の2点を示します。

トピック数（いくつのグループに分けるか）は予め決めておく
どのような意味を持つグループを作るかは予め人間が決めるのではなく、得られた確率から後で人間が解釈をする

特に後者についてはよく勘違いされている方がおりますのでご注意を。
例えば上記の例でいうと、トピック1は「サッカー」や「野球」という単語が高い確率で出現していることからトピック1は「スポーツ」と解釈を与えるといった感じです。
あくまでも、「統計学的に最適なグループ分けをした結果こうなりました」ということをお忘れなく。*1

2.2 LDAの理論

イメージだけでなく一応データサイエンティスト（見習い）感を出すために難しい話もします。

上記でイメージが掴めたなら読み飛ばしても大丈夫です。

では本題に戻り、早速論文のリンクを貼らせていただきます。

http://www.jmlr.org/papers/volume3/blei03a/blei03a.pdf

どこからどこまで話すべきなのか難しいですが、文書データは文書と単語の共起と捉えられ様々な生成モデルが研究されておりました。
この研究の一つとしてHoffmanはpLSAという確率的潜在モデル、いわゆる（確率的）トピックモデルが提案されました。
LDAはざっくり言うと、このpLSAを発展させ「潜在クラスの元での文書の出現確率（多項分布）に事前分布（ディリクレ分布）を仮定」したモデルとなっております。
相当端折って説明しているので、詳しくは論文のIntroductionを読んでください。

データの生成過程を表すグラフィカルモデルは以下のようになっております。
$P(\theta_d|\alpha)$ にディリクレ分布を仮定していることがミソです。 f:id:pira_nino:20180729154235p:plain

LDAのグラフィカルモデル（論文から引用）

このモデルの全文書（コーパス）の生成確率は以下の通りです。

f:id:pira_nino:20180729160246p:plain — 文書（コーパス）の生成確率（論文から引用）

この生成過程に対し、尤もらしいパラメータになるようにパラメータを推定します。*2

しかし、潜在クラスを仮定し、しかも事前分布を仮定しているために陽にパラメータ推定することができません。
そこで、Bleiらは「変分ベイズ」による推定を提案しました*3。

かなり駆け足で説明しましたが、LDAについての論文以外の日本語の説明としては以下の本に相当詳しく書かれております*4。

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

作者: 佐藤一誠,奥村学
出版社/メーカー: コロナ社
発売日: 2015/03/13
メディア: 単行本（ソフトカバー）
この商品を含むブログ (5件) を見る

3. 歌詞データへのLDAの適用

さて「文書データから単語のグルーピングを行う」LDAを歌詞データに適用します。

上記の数式を見るにめちゃくちゃ難しいモデルと感じますが、なんとgensimというパッケージで簡単に回せます。

radimrehurek.com

www.dskomei.com

3.1 コード

辞書の作成

使う単語のリスト化したものを辞書と呼び、gensimの型の辞書を作成します。

from gensim import corpora

text_list=list(df_all["英語の削除"])
dictionary =corpora.Dictionary(text_list)

実はこの段階でfilter_extremesを使うことで2文書(曲)以上に登場している単語のみを対象とするといったように前処理を行うことができます。

今回は

元のまま：7,763 語
1曲のみに：3,057語
2割以上の曲で出現する単語の削除後：3,028語
5割以上の曲で出現する単語の削除後：3,028語(2割と同じ)

という集計もあったので、「1曲のみに出現、2割以上の曲に出現する単語を削除」しました。*5

#1曲のみに出現、2割以上の曲に出現する単語を削除
dictionary.filter_extremes(no_below=2,no_above=0.2)
dictionary.save_as_text('最終辞書.txt')
# dictionary = corpora.Dictionary.load_from_text('最終辞書.txt')

ちなみにprint(dictionary.token2id)で辞書の中身を確認できます。

Corpusの作成

先ほど作成した辞書を元に、gensimのLDAの入力の型に文書データを変換しますよーってだけです。
ちょっと見慣れないdoc2bowという関数を使います。

corpus=[dictionary.doc2bow(tokens) for tokens in text_list]

LDAの学習

作成した辞書とコーパスを使ってLDAを回します。
といっても、gensimが便利すぎて下記のコードだけで学習ができます。

from gensim import models

#トピック数の設定
zk=20
#モデルの学習
model = gensim.models.LdaModel(corpus,
                               num_topics=zk,
                               id2word=dictionary,
                               random_state=2018
                              )
model.save('lda_bz.model')

これだけです。。。。gensimすごい。。。
ちなみに作者のmacbookでの334曲、一曲平均350単語くらいのコーパスの学習は1秒もかかりませんでした。

3.2 結果の確認

お待ちかねの結果です。
ここから先は（比較的）機械学習の知識がなくても読めます。

各トピックの出現確率

#各曲の各トピックへの所属確率の算出。(曲数×トピック数)のnumpy
Prob_songs=np.array(model.get_document_topics(corpus,minimum_probability=0))[:,:,1]

#DataFrameに収納
L=[ z for z in range(1,zk+1)]
col_name=list(map(lambda x: "Prob_"+str(x),L))
df_prob=pd.DataFrame(Prob_songs)
df_prob.columns=col_name

#所属確率最大のトピック番号も算出
df_prob["Max"]=df_prob.idxmax(axis=1)
def del_Prob(x):
    return int(x.split("_")[1])
df_prob["Max"]=df_prob["Max"].apply(lambda x : del_Prob(x))

#各トピックの出現確率の算出。曲のトピックへの所属確率を全曲で足して、全トピックで1になるように正規化
df_topic=pd.DataFrame(df_prob.drop("Max",axis=1).sum()/df_prob.drop("Max",axis=1).sum().sum())
df_topic.columns=["Prob"]
df_topic["Topic"]=[ z for z in range(1,zk+1)]

#可視化
plt.figure(figsize = (30,20))
ax= sns.barplot(x="Topic",y="Prob",data=df_topic,color="darkblue")
ax.set_xlabel("Topic",fontsize=50)
ax.set_ylabel("Prob",fontsize=50)
ax.tick_params(axis='x', labelsize=40)
ax.tick_params(axis='y', labelsize=40)
plt.show()

f:id:pira_nino:20180729183357p:plain — 各トピックの出現確率

特段1つのトピックの出現確率がぶっちぎっているなど、特に違和感のない結果になりました。*6

各トピックに高い確率で所属する単語

トピックごとの単語と出現確率のDataFrameを作成し、リストに格納
topic_word_prob=[]

for z in range(zk):
    word=[]
    prob=[]
    topic = model.show_topic(z,1000) #適当な単語数分

    for t in topic:
        word.append(t[0])
        prob.append(t[1])
        
    df_lda=pd.DataFrame({"word":word,"prob":prob})
    topic_word_prob.append(df_lda)

この結果を全トピック分表示するのも煩雑なので、数トピック分の所属確率の高いTop5の単語を以下の表に載せます。

Topic/Rank	1	2	3	4	5
Topic1	あなた	love	life	これ	世界
Topic7	baby	night	きみ	待っ	見つめ
Topic8	あなた	いつか	だれ	くる	like
Topic10	day	time	night	freedom	涙

率直な感想として「異なる意味の単語同士で綺麗に分かれた」結果が得られました。

これを元に独断と偏見で各トピックに解釈を与えます。

Topic 1	Topic 2	Topic 3	Topic 4
あなたと世界の果てでloveを欲しいlife	あなたと新しい世界でnight	旅立ちだれかと逢いひとつのloveの声聞いて	ラララあなたとDive Round

Topic 5	Topic 6	Topic 7	Topic 8
GUITAR KIDS RHAPSODY	知らないおまえとlove	nightにきみと見つめあってbaby	あなたといつか来るだれかのための声

Topic 9	Topic 10	Topic 11	Topic 12
baby yeah	Freedom Time day Now	ひとつaway	Lady Go Round

Topic 13	Topic 14	Topic 15	Topic 16
特に意味を持たない	あなたが欲しくてoh	特に意味を持たない	あなたと僕らでhey

Topic 17	Topic 18	Topic 19	Topic 20
今日あなたとbaby alright	特に意味を持たない	特に意味を持たない	特に意味を持たない

怒涛の20トピック分の解釈をしました。いい感じですね。

Word Cloudを用いた可視化

さて今回もWord Cloudを用いて可視化をします。
今回はトピックごとに所属確率の大きい単語を目立つように可視化することを目標とします。

データ入手編でも触れましたが、Word Cloudは単語の頻度の情報しか持ちません。そこで、トピックごとに単語の所属確率で重みづけた回数分の単語を生成させ、その擬似的な文書にWord Cloudを適用させます。

fig = plt.figure(figsize=(21, 12))

for z in range(zk):
    topic_text=[]
    for index, row in topic_word_prob[z].iloc[:2000,].iterrows():
        word=row["word"]
        weight=int(row["prob"]*1000)
        for _ in range(weight):
            topic_text.append(word)
            
    words = Counter(topic_text)
    bz_mask = np.array(Image.open("bz.png"))#適当な画像
    wc_bz = WordCloud(background_color="white",
                      max_words=3000,
                      max_font_size=70,
                      mask=bz_mask,
                      font_path=font_path
                     )
    wc_bz.generate_from_frequencies(words)

    plt.figure(figsize = (21,12))
    sns.set_style("whitegrid")
    plt.title('Topic_{}'.format(z+1), fontsize=30)
    plt.imshow(wc_bz, interpolation='bilinear')
    plt.axis("off")
    plt.savefig('20_Topic_{}.png'.format(z+1))
    plt.show()