ジョブレコメンデーションについてのリサーチまとめ

9月のイベントでジョブレコメンデーションについて調べて発表(実務と論文で学ぶジョブレコメンデーション最前線2022)しましたが、ブログの方が情報量が多いのと、最近のSlideShareが非常にみづらいものになっているのに加え、アップデートもしやすいのでこちらに随時残せる記事を残しておこうと思います。今日はクリスマスなので、誰かにとってはクリスマスプレゼントとなりうるでしょうか。

業界にいるものとしての思い

私は大学時代に経済学部に通っていたときに、労働市場の流動性に関して興味を持っていました。
スキルさえあればすぐに労働者が転職をするという競争的な労働市場について思い巡らすことがありました。

しかしながら、就職活動を通じて、取ってつけたような志望動機で大して勉強をしてこなかったであろう人々と同列に扱われ、「労働市場の効率性に関して、現実は不確実な情報のなかでエイヤで決めるしかないのかな」という印象を感じました。

就職活動では求人メディア系の会社に入社しました。スキルある人が高い給与を享受しやすい仕組み作りをデータ分析の力で実現できたらいいなと思って入りました。

ジョブマッチング界隈の課題

この業界の課題については推薦システム実践入門 ―仕事で使える導入ガイドの8章5節のドメインに応じた特徴と課題で、仕事の推薦について書かれているのも参考になりますので、そちらも合わせてお読みください。

双方の情報の充実化
- 情報の非対称性を払拭できるだけの情報がない中で企業も候補者も面接に臨んでいる
求人票の充実化
- 何かしらのエージェントを使う際に自分が本当に行きたい、行くべき職場がないかもしれない
求人票の質の向上
- 候補者が企業に直接応募するにしても、求人票に記されたレベル感が曖昧なため

ジョブマッチング界隈のデータの課題

スキルや職歴の構造化データ化
- 正確なスキル情報の記載
- 所属した企業に関しても法人番号などを付与できるといい
求人情報の記載項目の誤り是正、更新性アップ
- 今は求めていない要件があるとか、扱う技術に変化があったとかを反映
希望する仕事内容や職場の要望のデータ化
- 希望する仕事について構造化データになっていない

ジョブマッチング系の論文色々

・Job Recommender Systems: A Review(2021)
・Domain Adaptation for Resume Classification Using Convolutional Neural Networks(2017)
・Implicit Skills Extraction Using Document Embedding and Its Use in Job Recommendation(2020)
・Machine Learned Resume-Job Matching Solution(2016)
・SKILL: A System for Skill Identification and Normalization(2015)
・Effectiveness of job title based embeddings on résumé to job ad recommendation(2021)
・conSultantBERT: Fine-tuned Siamese Sentence-BERT for Matching Jobs and Job Seekers(2021)

Job Recommender Systems: A Review(2021)

論文の概要
- Job Recommender System(JRS)研究のサーベイ論文(2011〜2021)[計133本]
- 近年は求人などの大量のテキストを用いたジョブレコメンデーションが多い。深層学習も使われている。BERTもCNNも使われている。ハイブリッドモデルが多い傾向。
- 近年ではアルゴリズムの公平性(年齢/性別での差別)に関心が持たれており、特徴量として差別的なものがなくても学習してしまう可能性が指摘されている。
- 求職者とリクルーターのインタラクティブなデータ(click/skip)の利活用が重要となっている。研究者はRecsysなどのコンペのデータに依存しており、企業がデータセットの提供をするのが限定的であるのが問題としてある。JRSの発展のためにはKaggleやRecsysのコンペが重要とされている。
参考になるところ
- 求職者と採用担当者のインタラクティブなデータはこのコミュニティでデータとして貴重であること
- 最近はTransformer系のアプローチが使われ始めていること

Domain Adaptation for Resume Classification Using Convolutional Neural Networks(2017)

論文の概要
- 27の職種のラベルが付与された求人票のデータ（テキスト：短文）をもとにCNNで分類器を作成
- ラベルがない職務経歴データに関しては求人票で学習したモデルで推論
- fastTextでの分類よりも性能が向上した
- 求人票の数 > 越えれない壁 > 職務経歴データの数
- CNNで学習する際の特徴量は学習済みのWord2Vecによる単語の分散表現で、同じ長さのものを使っている。
- 最大値プーリングでの重要な特徴抽出、ドロップアウトの正則化など画像分類タスクと同様になされ、soft-max関数でマルチクラスの確率を出力している。
- 苦手な職種でのファインチューニングがfuture work
参考になるところ
- 職務経歴書と違って収集が容易な求人票のデータをもとに学習し、職務経歴に関するラベルを推論をするスタイル
- テキストデータに対してCNNを使っている

Implicit Skills Extraction Using Document Embedding and Its Use in Job Recommendation(2020)

論文の概要
- 職務経歴書や求人からNLPでスキル表現を抽出
- 暗黙的なスキル表現も抽出
- 候補者と求人のマッチングのスコア(Affinity Score)を定義
- 暗黙的スキル情報を使うとマッチングの性能がアップすることがわかった
- テキストに対して3つの観点でスキル抽出
  ・NER(Named Entity Recognition(固有表現抽出))を適用→文書内のスキルの位置を見つける、カテゴライズを行う
  ・PoS(Part of Speech(品詞の付与))を適用→スキル表現かどうかは専門家によるアノテーションを実施
  ・色々集めた辞書にあるか計算→Wikipedia、Onet/Hopeというサイトなどで自動での拡張型スキル辞書用意
- 暗黙的なスキル表現も抽出
  暗黙スキル（Web DevelopmentやHard Workingなどの表現）の考慮→分散表現の平均値で表現
- スキルとの関連度スコアの定義
  職務経歴書や求人からのスキル抽出のためにスコアを用いる
- 表現の除外
  関連度スコアに閾値を設けて下回ったらその表現を除外するなどを行う。
- 関連度スコア(x)
  ・NERで抽出した表現の信頼度(S1)
  ・文法が、事前に定義したルールにどれだけ合っているか(S2)
  ・そのスキルがNERで得たスキルっぽいもの、ポテンシャルスキル、辞書スキルの3つの辞書にどれだけ合致したかの割合(S3)
  ・その単語とスキル辞書との分散表現における最大のcos類似度の値(S4)
  ・パラメータは経験則的に決める
  ・関連スコアは0~1の間をとる
- アフィニティスコア
  ・候補者と求人のマッチの適切さの指標
  ・エッジウェイトの平均としている
- エッジウェイト(Y)
  ・スキル間のcos類似度(E1)
  ・文書を通じたスキルの出現頻度の割合(E2)
  ・明示的なスキルがあるかどうか(E3)→暗黙だと減点
  ・ウェイト(ω)はそのマッチングが成功したかどうかで重み付ける。
  ・関連スコアは0~1の間をとる

参考になるところ

暗黙的なスキルに関して着目したところ
スコアリングのモデルを定め、手元にあるデータとWikipediaなどの集合知も使ってかなり慎重にスキル抽出を行っている

Machine Learned Resume-Job Matching Solution(2016)

論文の概要
- 深層学習(CNN)を用いて、リッチな特徴量を用意した上で、職務経歴書と求人のマッチング最適化を行っている
- 職務経歴書の数は4.7万枚
- KerasやScikit-learnで書かれており軽量
- 浅いモデルとしてXGBoostを、深いモデルとしてCNNの両方を、最終的に多数決によるアンサンブルにすることで性能が向上
- アンサンブル手法としてはバギングを行っている
- リッチな特徴量
  ・Manual Feature(年齢とか性別とか)
  ・Cluster Feasture(クラスタリングしたもの、トピック(LDA))
  ・Semantic Feature(意味的な類似度)
- モデルはGPUでもCPUでも動く
- ソースコードはGitHubで公開されてる
参考になるところ
- CNNながらGPUなしでも計算できるようなモデルのアーキテクチャになっている
- 特徴量はマニュアルで作成したものと自動で抽出したものの両方をうまく使う

SKILL: A System for Skill Identification and Normalization(2015)

論文の概要
- 職務経歴書や求人に記されているスキルの抽出に関するNER(Named Entity Recognition：固有表現抽出)の研究
- 老舗の求人サイトであるCareerBuilderの
  社員による研究
- Wikipediaベースのカテゴリ分類、Google検索による集合知を利用してスキル表現の抽出をしている
- Word2Vecを用いて、あるスキルと関連しているスキルをもとにタグ付をするアプローチが紹介されている
- データ
  大量の職務経歴書と求人票
- アプローチ
  ・スキル分類→wikipediaベースのAPI、国の定めた定義など
  ・スキル判定→単一語(そのまま利用)、複数語(Google検索結果から判断(集合知))
  ・スキルタグ付け→Word2Vecをユニグラムで学習し、関連度を計算しスキルのタグ付け
参考になるところ
- Google検索を使って得た集合知をスキルの分類に活かしている
- Word2vecで計算した類似度をもとにスキルのタグ付けをしている

Effectiveness of job title based embeddings on résumé to job ad recommendation(2021)

論文の概要
- 求人票の表現に対して、分散表現を用いた求人レコメンドの研究
- 求人票全体を使うよりも、タイトル部分を使ったほうがレコメンドの性能が高くなった
  ・分散表現はWord2VecやDoc2Vecの学習済みや学習したものを色々試している
  ・Cold-Start問題に対しても分散表現は有効な手段であるとしている
参考になるところ
- 汚い本文のテキストを使うよりもタイトルを使ったほうが良いと言う、岡目八目というか、労せずに成果を出す作法を知れたという気持ち
- 分散表現がCold-Start問題に対する有効な手段という可能性

conSultantBERT: Fine-tuned Siamese Sentence-BERT for Matching Jobs and Job Seekers(2021)

論文の概要
- データは27万件で正例が12.6万件、負例が10.9万件で職務経歴書（ユニークだと15万件）と求人票（2.3万件）からなる。ラベルはキャリアアドバイザーが付けたもの。
- 各センテンスに対してBERTを適用し、文書単位で平均を取るなどしている。文書は最初の512トークンを使うのがよいとしている。
- BERTを分類問題・回帰問題についてそれぞれSiamese network（シャムネットワーク）を用いてファインチューニング
  ・出力層は求職者のレジュメと求人票のCOS類似度をマッチングスコアとしている
  ・埋め込み層のデータで教師あり学習をランダムフォレストで行い、COS類似度ではなくランダムフォレストの吐いたスコアを使って推論するバージョンも用意している
- 回帰タスクでのファインチューニングが一番よかった
参考になるところ
- BERTをセンテンスに対して使って、平均値をとりそれを特徴量にしている
- 全文を使わないほうが分類性能が良い可能性
- 埋め込み層のデータで教師あり学習をして、それのスコアを用いるというアプローチ

BLSTMを用いた文書分類でデザイナーズマンション予測に再挑戦

はじめに

仕事で深層学習を使うことはないのですが、扱える技術の幅を広げていきたいなと思い、BLSTMを用いた文書分類についてkerasでの簡単なサンプルコードをやってみようと思います。データは以前集めた、某不動産紹介サイトの賃貸マンションの設備に関する文書とそれがデザイナーズマンションかどうかのラベルを使います。そして文書内の単語からその文書がデザイナーズマンションかどうかを予測します。前回はAUCで83%だったので、それを超えれると良いですね。

単純なRNNとは

モチベーション
- フィードフォワード型のニューラルネットワークではうまく扱うことができない時系列データをうまく扱えるようにすること。
特徴
- 入力が互いに関係している（多層パーセプトロンの際は置かれていない仮定）
  - 直訳すると循環するニューラルネットワークとなる。
- 最初の文の単語が2つ目、3つ目の単語に影響を与える可能性を考慮。

具体的な関数の形としては、 $$ h_t = \tanh (W h_{t-1} + U x_t) \\\ y_t = softmax(V h_t) $$ で与えられる。

$ h_t $は隠れ層の状態を、$ x_t $は入力変数を、$ y_t $は出力ベクトルを表している。（他のアルファベットは重み行列）

LSTMとは

モチベーション
- 単純なRNNの勾配消失問題を解決するために提案された手法。
特徴
- 単純なRNNに置き換えることで性能が大幅に向上することも珍しくない。
  - 時系列データ、長い文章、録音データからなる長期的なパターンを取り出すことを得意としている手法。
- 勾配消失問題に強い
- Long Short-Term Memory：長短期記憶
  - 長期依存性を学習できるRNNの亜種。
  - 入力ゲート、忘却ゲート、出力ゲート、内部隠れ層という4つの層が相互に関わり合う。重要な入力を認識し、それを長期状態に格納すること、必要な限りで記憶を保持すること、必要なときに記憶を取り出すことを学習する。
    - 入力ゲート：内部隠れ層のどの部分を長期状態に加えるかを決める。
    - 忘却ゲート：長期状態のどの部分を消去するか決める。
    - 出力ゲート：各タイムステップで、長期状態のどの部分を読み出し、出力するかを決める。

$$
i = \sigma (W_i h_{t-1} + U_i x_t) \\
f = \sigma (W_f h_{t-1} + U_f x_t) \\
o = \sigma (W_o h_{t-1} + U_ox_t ) \\
g = \tanh (W_g h_{t-1} + U_g x_t) \\
c_t = (c_{t-1} \otimes f ) \otimes ( g \otimes i) \\
h_t = \tanh (c_t) \otimes o
$$ i：入力ゲート(input)
f：忘却ゲート(forget)
o：出力ゲート(output)
$ \sigma $：シグモイド関数
g：内部隠れ層状態
$c_t $：時刻tにおけるセル状態
$h_t $：時刻tにおける隠れ状態

Bidirectional LSTMとは

モチベーション
- 従来のRNNの制約を緩和するために導入。
特徴
- ある特定の時点で過去と将来の利用可能な入力情報を用いて訓練するネットワーク
  - 従来のRNNのニューロンをフォワード(未来)なものとバックワード(過去)なものとの2つに分ける。
    - 2つのLSTMを訓練している。
  - 全体のコンテキストを利用して推定することができるのが良いらしい。
    - 文章で言うと、文章の前(過去)と後(未来)を考慮して分類などを行うイメージ。
  - 人間もコンテキストを先読みしながら解釈することもある。
- BLSTMは全てのシークエンスの予測問題において有効ではないが、適切に扱えるドメインで良い結果を残す。
- 1997年とかなり歴史があるものらしい

（出典：Deep Dive into Bidirectional LSTM）

Bidirectional LSTMで文書分類

今回は、BLSTMを使って、マンションの設備に関するテキスト情報から、そのマンションがデザイナーズマンションかどうかを予測します。全体のソースコードはGoogle Colabを御覧ください。

データ

import pandas as pd
from sklearn.model_selection import train_test_split

# GitHubにおいてあるマンションのデータを読み込む
raw_data = pd.read_csv("https://github.com/KamonohashiPerry/kamonohashiperry.com/blob/master/designers_apartment/designers_apartment.csv?raw=true")
raw_data["text"] = raw_data["text"].str.replace("デザイナーズ", "") 

# 訓練データとテストデータを分割
train_x, test_x, train_y, test_y = train_test_split(raw_data["text"], raw_data["designer_flag"], test_size=0.3, random_state=123)

import pandas as pd

from sklearn.model_selection import train_test_split

# GitHubにおいてあるマンションのデータを読み込む

raw_data = pd.read_csv("https://github.com/KamonohashiPerry/kamonohashiperry.com/blob/master/designers_apartment/designers_apartment.csv?raw=true")

raw_data["text"] = raw_data["text"].str.replace("デザイナーズ", "")

# 訓練データとテストデータを分割

train_x, test_x, train_y, test_y = train_test_split(raw_data["text"], raw_data["designer_flag"], test_size=0.3, random_state=123)

データを確認してみると、

train_x
246     バストイレ別 バルコニー エアコン クロゼット フローリング 浴室乾燥機 オートロック シス...
1402    バストイレ別 バルコニー エアコン ガスコンロ対応 クロゼット フローリング TVインターホ...
1586    バストイレ別 バルコニー エアコン ガスコンロ対応 クロゼット フローリング TVインターホ...
375     バストイレ別 バルコニー エアコン クロゼット フローリング 浴室乾燥機 オートロック 室内...
1272    エアコン フローリング 室内洗濯置 シューズボックス 角住戸 光ファイバー 2面採光 最上階...
1122    バストイレ別 バルコニー エアコン ガスコンロ対応 クロゼット フローリング TVインターホ...
1346    バストイレ別 バルコニー エアコン ガスコンロ対応 クロゼット フローリング オートロック...
1406    バストイレ別 バルコニー エアコン ガスコンロ対応 フローリング シャワー付洗面台 オートロ...
1389    バストイレ別 エアコン クロゼット フローリング シューズボックス 押入 光ファイバー 即入...
1534    バルコニー エアコン クロゼット TVインターホン オートロック 陽当り良好 シューズボック...

train_x

246 バストイレ別バルコニーエアコンクロゼットフローリング浴室乾燥機オートロックシス...

1402 バストイレ別バルコニーエアコンガスコンロ対応クロゼットフローリング TVインターホ...

1586 バストイレ別バルコニーエアコンガスコンロ対応クロゼットフローリング TVインターホ...

375 バストイレ別バルコニーエアコンクロゼットフローリング浴室乾燥機オートロック室内...

1272 エアコンフローリング室内洗濯置シューズボックス角住戸光ファイバー 2面採光最上階...

1122 バストイレ別バルコニーエアコンガスコンロ対応クロゼットフローリング TVインターホ...

1346 バストイレ別バルコニーエアコンガスコンロ対応クロゼットフローリングオートロック...

1406 バストイレ別バルコニーエアコンガスコンロ対応フローリングシャワー付洗面台オートロ...

1389 バストイレ別エアコンクロゼットフローリングシューズボックス押入光ファイバー即入...

1534 バルコニーエアコンクロゼット TVインターホンオートロック陽当り良好シューズボック...

train_y

246 1

1402 0

1586 0

375 1

1272 0

1122 0

1346 0

1406 0

1389 0

1534 0

今回のデータがテキストとラベルからなっていることがわかる。なお、データ数は1864件あり、そのうち16%がデザイナーズマンションです。

前処理

Google Colab上で形態素解析を行うために、MeCabをインストールする。

!apt-get -q -y install sudo file mecab libmecab-dev mecab-ipadic-utf8 git curl python-mecab
# せっかくなので、Neologdを使えるようにする。
!git clone --depth 1 https://github.com/neologd/mecab-ipadic-neologd.git
!echo yes | mecab-ipadic-neologd/bin/install-mecab-ipadic-neologd -n
!cat /etc/mecabrc
!sed -e "s!/var/lib/mecab/dic/debian!/usr/lib/x86_64-linux-gnu/mecab/dic/mecab-ipadic-neologd!g" /etc/mecabrc > /etc/mecabrc.new
!cat /etc/mecabrc.new
!cp /etc/mecabrc /etc/mecabrc.org
!cp /etc/mecabrc.new /etc/mecabrc
!cat /etc/mecabrc
!apt-get -q -y install swig
!pip install mecab-python3

!apt-get -q -y install sudo file mecab libmecab-dev mecab-ipadic-utf8 git curl python-mecab

# せっかくなので、Neologdを使えるようにする。

!git clone --depth 1 https://github.com/neologd/mecab-ipadic-neologd.git

!echo yes | mecab-ipadic-neologd/bin/install-mecab-ipadic-neologd -n

!cat /etc/mecabrc

!sed -e "s!/var/lib/mecab/dic/debian!/usr/lib/x86_64-linux-gnu/mecab/dic/mecab-ipadic-neologd!g" /etc/mecabrc > /etc/mecabrc.new

!cat /etc/mecabrc.new

!cp /etc/mecabrc /etc/mecabrc.org

!cp /etc/mecabrc.new /etc/mecabrc

!cat /etc/mecabrc

!apt-get -q -y install swig

!pip install mecab-python3

これでMaCab NeologdをGoogle ColabのPythonで実行できます。

テキストデータの前処理を行うために名詞だけを抽出してリスト形式で返す関数を定義します。

# 形態素解析を行い、名詞だけ抽出し、かつ名詞の中でも活用を絞り込んで抽出する関数
def nouns_extract(line):
    import re
    import MeCab
    keyword=[]
    m = MeCab.Tagger(' -d  /usr/lib/x86_64-linux-gnu/mecab/dic/mecab-ipadic-neologd ')
    for l in m.parse(line).splitlines():
        if l != 'EOS' and l.split('\t')[1].split(',')[0] == '名詞':
            if l != 'EOS' and re.search('^一般$|^形容動詞語幹$|^サ変接続$|^副詞可能$|^固有名詞

# 形態素解析を行い、名詞だけ抽出し、かつ名詞の中でも活用を絞り込んで抽出する関数

def nouns_extract(line):

import re

import MeCab

keyword=[]

m = MeCab.Tagger(' -d /usr/lib/x86_64-linux-gnu/mecab/dic/mecab-ipadic-neologd ')

for l in m.parse(line).splitlines():

if l != 'EOS' and l.split('\t')[1].split(',')[0] == '名詞':

if l != 'EOS' and re.search('^一般$|^形容動詞語幹$|^サ変接続$|^副詞可能$|^固有名詞

ここで、語彙に関するデータを作成します。

import collections

# 語彙に関するデータを作成するための処理
maxlen = 0
word_freqs = collections.Counter()
num_recs = 0

for sentence in list(train_x):
  words = nouns_extract(sentence)
  maxlen = max(maxlen, len(words)) # リストの最大値を取得
  # 単語の出現頻度を集計
  for word in words:
    word_freqs[word] += 1
  num_recs += 1

import collections

# 語彙に関するデータを作成するための処理

maxlen = 0

word_freqs = collections.Counter()

num_recs = 0

for sentence in list(train_x):

words = nouns_extract(sentence)

maxlen = max(maxlen, len(words)) # リストの最大値を取得

# 単語の出現頻度を集計

for word in words:

word_freqs[word] += 1

num_recs += 1

続いて、単語とそのインデックスからなるdict形式のデータを作成します。

# 語彙のなかから2000語を選ぶ
MAX_FEATURES = 2000
# 文の長さの最大値
MAX_SENTENCE_LENGTH = 40
# 語彙数
vocab_size = min(MAX_FEATURES, len(word_freqs)) + 2
# 単語とそのインデックスからなるdict形式のデータ
word2index = {x[0]: i+2 for i, x in enumerate(word_freqs.most_common(MAX_FEATURES))}

# 0番目と1番目にパディング用(PAD)と擬似単語(UNK)を指定する
# パディングとは系列の長さを揃えること。
word2index['PAD'] = 0
word2index['UNK'] = 1

# インデックスとその単語からなるdict形式のデータ
index2word = {v:k for k, v in word2index.items()}

# 語彙のなかから2000語を選ぶ

MAX_FEATURES = 2000

# 文の長さの最大値

MAX_SENTENCE_LENGTH = 40

# 語彙数

vocab_size = min(MAX_FEATURES, len(word_freqs)) + 2

# 単語とそのインデックスからなるdict形式のデータ

word2index = {x[0]: i+2 for i, x in enumerate(word_freqs.most_common(MAX_FEATURES))}

# 0番目と1番目にパディング用(PAD)と擬似単語(UNK)を指定する

# パディングとは系列の長さを揃えること。

word2index['PAD'] = 0

word2index['UNK'] = 1

# インデックスとその単語からなるdict形式のデータ

index2word = {v:k for k, v in word2index.items()}

最後に、これまでに作成した語彙とそのインデックスのデータを用いて、実際に学習に使う入力データと出力データを作成します。

import numpy as np

# 空の変数を作成
X = np.empty((num_recs, ), dtype=list)
y = np.zeros((num_recs, ))
i = 0

train_x_list = list(train_x)
train_y_list = list(train_y)

for i in range(len(train_x_list)):
  label = train_y_list[i]
  words = nouns_extract(train_x_list[i])
  # 特徴量用の空の変数
  seqs = []
  for word in words:
    # 単語を照合してインデックスを取得
    if word in word2index:
      seqs.append(word2index[word])
    else:
      # 語彙に存在しない単語の場合はUNKのインデックスを付与する          
      seqs.append(word2index["UNK"])

  X[i] = seqs
  y[i] = int(label)
  i += 1

# パディングして、同じ長さに揃えてくれる関数
X = sequence.pad_sequences(X, maxlen=MAX_SENTENCE_LENGTH)

import numpy as np

# 空の変数を作成

X = np.empty((num_recs, ), dtype=list)

y = np.zeros((num_recs, ))

i = 0

train_x_list = list(train_x)

train_y_list = list(train_y)

for i in range(len(train_x_list)):

label = train_y_list[i]

words = nouns_extract(train_x_list[i])

# 特徴量用の空の変数

seqs = []

for word in words:

# 単語を照合してインデックスを取得

if word in word2index:

seqs.append(word2index[word])

else:

# 語彙に存在しない単語の場合はUNKのインデックスを付与する

seqs.append(word2index["UNK"])

X[i] = seqs

y[i] = int(label)

i += 1

# パディングして、同じ長さに揃えてくれる関数

X = sequence.pad_sequences(X, maxlen=MAX_SENTENCE_LENGTH)

実行

先程作成したデータを訓練データとテストデータに分けます。

Xtrain, Xtest, ytrain, ytest = train_test_split(X, y, test_size=0.2, random_state=42)

1	Xtrain, Xtest, ytrain, ytest = train_test_split(X, y, test_size=0.2, random_state=42)

ここで、AUCを計算するための関数を定義します。（まるまる拝借しました）

import tensorflow as tf
from sklearn.metrics import roc_auc_score
from keras.callbacks import Callback, EarlyStopping

# define roc_callback, inspired by https://github.com/keras-team/keras/issues/6050#issuecomment-329996505
def auc_roc(y_true, y_pred):
    # any tensorflow metric
    value, update_op = tf.contrib.metrics.streaming_auc(y_pred, y_true)
    
    # find all variables created for this metric
    metric_vars = [i for i in tf.local_variables() if 'auc_roc' in i.name.split('/')[1]]
    
    # Add metric variables to GLOBAL_VARIABLES collection.
    # They will be initialized for new session.
    for v in metric_vars:
        tf.add_to_collection(tf.GraphKeys.GLOBAL_VARIABLES, v)

    # force to update metric values
    with tf.control_dependencies([update_op]):
        value = tf.identity(value)
        return value

import tensorflow as tf

from sklearn.metrics import roc_auc_score

from keras.callbacks import Callback, EarlyStopping

# define roc_callback, inspired by https://github.com/keras-team/keras/issues/6050#issuecomment-329996505

def auc_roc(y_true, y_pred):

# any tensorflow metric

value, update_op = tf.contrib.metrics.streaming_auc(y_pred, y_true)

# find all variables created for this metric

metric_vars = [i for i in tf.local_variables() if 'auc_roc' in i.name.split('/')[1]]

# Add metric variables to GLOBAL_VARIABLES collection.

# They will be initialized for new session.

for v in metric_vars:

tf.add_to_collection(tf.GraphKeys.GLOBAL_VARIABLES, v)

# force to update metric values

with tf.control_dependencies([update_op]):

value = tf.identity(value)

return value

続いて、kerasを用いて、ネットワークを構築し、コンパイルを行います。

from keras.layers import Activation, Dense, Dropout, Embedding, LSTM, Bidirectional
from keras.models import Sequential
from keras.preprocessing import sequence
from keras.callbacks import TensorBoard

LOG_DIR = "drive/Colab Notebooks/NLP/logs"

# 分散表現層の次元
EMBEDDING_SIZE = 128
# 隠れ層の数
HIDDEN_LAYER_SIZE = 64
# バッチサイズ
BATCH_SIZE = 32
# エポックの数
NUM_EPOCHS = 10

model = Sequential()

# 分散表現の層
model.add(Embedding(vocab_size, EMBEDDING_SIZE, input_length=MAX_SENTENCE_LENGTH))
# Bidirectional LSTMの層
model.add(Bidirectional(LSTM(HIDDEN_LAYER_SIZE )))
# ドロップアウト層
model.add(Dropout(0.5))
# 活性化層
model.add(Dense(1, activation='sigmoid'))

# モデルのコンパイル
model.compile(loss="binary_crossentropy", optimizer="adam", metrics=["accuracy", auc_roc])

my_callbacks = [EarlyStopping(monitor='auc_roc', patience=300, verbose=1, mode='max'), TensorBoard(LOG_DIR)]

# ネットワークの学習
history = model.fit(Xtrain, ytrain, batch_size=BATCH_SIZE,
                    epochs=NUM_EPOCHS,
                    callbacks=my_callbacks,
                    validation_data=(Xtest, ytest))

from keras.layers import Activation, Dense, Dropout, Embedding, LSTM, Bidirectional

from keras.models import Sequential

from keras.preprocessing import sequence

from keras.callbacks import TensorBoard

LOG_DIR = "drive/Colab Notebooks/NLP/logs"

# 分散表現層の次元

EMBEDDING_SIZE = 128

# 隠れ層の数

HIDDEN_LAYER_SIZE = 64

# バッチサイズ

BATCH_SIZE = 32

# エポックの数

NUM_EPOCHS = 10

model = Sequential()

# 分散表現の層

model.add(Embedding(vocab_size, EMBEDDING_SIZE, input_length=MAX_SENTENCE_LENGTH))

# Bidirectional LSTMの層

model.add(Bidirectional(LSTM(HIDDEN_LAYER_SIZE )))

# ドロップアウト層

model.add(Dropout(0.5))

# 活性化層

model.add(Dense(1, activation='sigmoid'))

# モデルのコンパイル

model.compile(loss="binary_crossentropy", optimizer="adam", metrics=["accuracy", auc_roc])

my_callbacks = [EarlyStopping(monitor='auc_roc', patience=300, verbose=1, mode='max'), TensorBoard(LOG_DIR)]

# ネットワークの学習

history = model.fit(Xtrain, ytrain, batch_size=BATCH_SIZE,

epochs=NUM_EPOCHS,

callbacks=my_callbacks,

validation_data=(Xtest, ytest))

テストデータでの予測精度を確認します。

# testデータでの予測
score, acc, auc = model.evaluate(Xtest, ytest, batch_size=BATCH_SIZE)
print("Test score: {:.3f}, accuracy: {:.3f}, auc: {:.3f}".format(score, acc, auc))

Test score: 0.444, accuracy: 0.824, auc: 0.844

# testデータでの予測

score, acc, auc = model.evaluate(Xtest, ytest, batch_size=BATCH_SIZE)

print("Test score: {:.3f}, accuracy: {:.3f}, auc: {:.3f}".format(score, acc, auc))

Test score: 0.444, accuracy: 0.824, auc: 0.844

実際のラベルとテキストデータを見比べてみます。

for i in range(20):
  idx = np.random.randint(len(Xtest))
  xtest = Xtest[idx].reshape(1, 40)
  ylabel = ytest[idx]
  ypred = model.predict(xtest)[0][0]
  sent = " ".join([index2word[x] for x in xtest[0].tolist() if x != 0 ])
  print("{:.0f}\t{:.0f}\t{}".format(ypred, ylabel, sent))

0   0   バルコニー エアコン クロゼット フローリング 室内 洗濯 置 シューズボックス システムキッチン 住戸 外壁 タイル 入居 ガスレンジ 都市ガス
0   0   ウォークインクロゼット 保証 不要 敷金 バイク 置場 ネット 専用 回線 トランクルーム 事務所 相談 楽器 相談 トイレ ルームシェア 相談 利用 沿線 利用 徒歩 10分 専有 面積 平面 駐車場 LDK 居室 都市ガス 高速 ネット 対応 礼金 保証 会社 利用 初期 費用 カード 決済
0   0   バス トイレ エアコン ガスコンロ 対応 クロゼット フローリング 室内 洗濯 置 シューズボックス システムキッチン 脱衣所 入居 礼金 不要 最上階 敷金 ガスレンジ ダブル ロック キー 上 都市ガス
0   0   室内 洗濯 置 焚 機能 浴室 住戸 脱衣所 洗面所 独立 駐輪場 光ファイバー 入居 ペット 相談 敷金 居室 フローリング 採光 収納 内装 リフォーム 済 事務所 相談 ルームシェア 相談 利用 沿線 利用 徒歩 5分 東南 都市ガス 玄関 収納 礼金 保証 会社 利用
0   0   居室 フローリング 採光 沿線 利用 内装 リフォーム 済 眺望 良好 事務所 相談 24時間 換気 システム 南面 リビング 外装 リフォーム 済 利用 利用 沿線 利用 徒歩 10分 敷地内 ごみ 置き場 都市ガス BS 保証 会社 利用 初期 費用 カード 決済 通風 良好
0   0   TV インターホン 浴室 乾燥機 オートロック 室内 洗濯 置 システムキッチン 住戸 エレベーター 洗面所 独立 コンロ 宅配ボックス 光ファイバー 外壁 タイル BS CS 防犯カメラ 保証 不要 敷金 保証金 不要 プール 利用 沿線 利用 徒歩 10分 24時間 ゴミ 敷地内 ごみ 置き場 5年 都市ガス 礼金
0   0   バルコニー エアコン クロゼット フローリング 室内 洗濯 置 シューズボックス システムキッチン 住戸 外壁 タイル 入居 ガスレンジ 都市ガス
0   1   バス トイレ バルコニー エアコン クロゼット フローリング TV インターホン 浴室 乾燥機 オートロック 室内 洗濯 置 シューズボックス システムキッチン 温水洗浄便座 エレベーター 洗面所 独立 コンロ 駐輪場 礼金 不要 バイク 置場 ガスレンジ ネット 専用 回線 利用 沿線 利用 徒歩 10分 都市ガス BS
1   1   システム 居室 フローリング ディンプルキー ダブル ロック キー 24時間 換気 システム フロント サービス 耐火構造 耐震 構造 利用 沿線 利用 徒歩 5分 24時間 ゴミ 敷地内 ごみ 置き場 セキュリティ 会社 加入 済 間接照明 LAN 保証 会社 利用 初期 費用 カード 決済 通風 良好
0   0   バス トイレ バルコニー フローリング TV インターホン 浴室 乾燥機 オートロック 室内 洗濯 置 システムキッチン 住戸 エレベーター 洗面所 独立 CATV 入居 防犯カメラ グリル 居室 洋室 仲 エアコン 納戸 徒歩 10分 専有 面積 南西
0   0   エアコン ガスコンロ 対応 クロゼット フローリング シャワー 洗面 TV インターホン 浴室 乾燥機 室内 洗濯 置 シューズボックス 温水洗浄便座 脱衣所 洗面所 独立 駐輪場 CATV 保証 不要 CATV インターネット 平坦 保証 不要 始発駅 利用 沿線 利用 徒歩 10分 バス停 徒歩 3分 南西 都市ガス BS
0   1   住戸 浴室 窓 24時間 換気 システム クロゼット ルームシェア 相談 南面 リビング 学生 相談 一部 フローリング 利用 沿線 利用 徒歩 5分 駅 徒歩 10分 24時間 ゴミ 敷地内 ごみ 置き場 都市ガス 洗面所 ドア 南面 バルコニー BS 初期 費用 カード 決済 通風 良好0 0   対面 キッチン IH クッキング ヒーター グリル ウォークインクロゼット 保証 不要 敷金 二人 入居 相談 ネット 使用 不要 テラス 保証 不要 3年 南面 リビング 浴室 利用 沿線 利用 バス停 徒歩 3分 敷地内 ごみ 置き場 LDK 都市ガス 南面 バルコニー 天井 シューズ クロゼット 礼金
0   0   バス トイレ バルコニー エアコン フローリング 室内 洗濯 置 シューズボックス 焚 機能 浴室 駐輪場 宅配ボックス 入居 礼金 不要 敷金 不要 IH クッキング ヒーター エレベーター 2基 利用 沿線 利用 徒歩 5分 都市ガス
1   1   コンロ 宅配ボックス 光ファイバー 外壁 タイル BS CS 防犯カメラ 照明 保証人 不要 敷金 沿線 利用 ディンプルキー 洗濯機 耐火構造 始発駅 沿線 利用 徒歩 10分 24時間 ゴミ 敷地内 ごみ 置き場 居室 都市ガス 高速 ネット 対応 礼金 保証 会社 利用 初期 費用 カード 決済
0   0   バルコニー エアコン クロゼット フローリング シューズボックス 住戸 エレベーター 入居 防犯カメラ 分譲 賃貸 仲 眺望 良好 保証 不要 利用 沿線 利用 徒歩 5分 敷地内 ごみ 置き場 日勤 管理 敷金 礼金 不要 通風 良好
1   1   24時間 換気 システム 南面 リビング 一部 フローリング ISDN 対応 耐火構造 耐震 構造 利用 徒歩 5分 バス停 徒歩 3分 制震 構造 高層 24時間 ゴミ 敷地内 ごみ 置き場 都市ガス 南面 バルコニー 玄関 収納 BS 高速 ネット 対応 LAN 初期 費用 カード 決済
0   0   バイク 置場 南西 住戸 ガスレンジ ディンプルキー 内装 リフォーム 済 眺望 良好 照明 センサー 利用 沿線 利用 徒歩 10分 24時間 ゴミ 敷地内 ごみ 置き場 南西 当社 管理 物件 都市ガス 洗面所 ドア 南面 バルコニー 室内 物干 礼金 保証 会社 利用 通風 良好
0   0   駅 平坦 ネット 使用 不要 ダブル ロック キー 保証 不要 2年 24時間 換気 システム 複層ガラス 照明 センサー 利用 徒歩 10分 外壁 サイディング 24時間 ゴミ 敷地内 ごみ 置き場 南西 都市ガス 玄関 収納 年内 入居 年度内 入居 礼金 初期 費用 カード 決済
0   1   フローリング TV インターホン 浴室 乾燥機 オートロック 室内 洗濯 置 陽 良好 シューズボックス システムキッチン 焚 機能 浴室 住戸 エレベーター 洗面化粧台 コンロ 駐輪場 宅配ボックス CATV 光ファイバー BS CS 対面 キッチン 沿線 利用 楽器 相談 ウッドデッキ 利用 徒歩 10分 東南 LDK 居室 都市ガス

for i in range(20):

idx = np.random.randint(len(Xtest))

xtest = Xtest[idx].reshape(1, 40)

ylabel = ytest[idx]

ypred = model.predict(xtest)[0][0]

sent = " ".join([index2word[x] for x in xtest[0].tolist() if x != 0 ])

print("{:.0f}\t{:.0f}\t{}".format(ypred, ylabel, sent))

0 0 バルコニーエアコンクロゼットフローリング室内洗濯置シューズボックスシステムキッチン住戸外壁タイル入居ガスレンジ都市ガス

0 0 ウォークインクロゼット保証不要敷金バイク置場ネット専用回線トランクルーム事務所相談楽器相談トイレルームシェア相談利用沿線利用徒歩 10分専有面積平面駐車場 LDK 居室都市ガス高速ネット対応礼金保証会社利用初期費用カード決済

0 0 バストイレエアコンガスコンロ対応クロゼットフローリング室内洗濯置シューズボックスシステムキッチン脱衣所入居礼金不要最上階敷金ガスレンジダブルロックキー上都市ガス

0 0 室内洗濯置焚機能浴室住戸脱衣所洗面所独立駐輪場光ファイバー入居ペット相談敷金居室フローリング採光収納内装リフォーム済事務所相談ルームシェア相談利用沿線利用徒歩 5分東南都市ガス玄関収納礼金保証会社利用

0 0 居室フローリング採光沿線利用内装リフォーム済眺望良好事務所相談 24時間換気システム南面リビング外装リフォーム済利用利用沿線利用徒歩 10分敷地内ごみ置き場都市ガス BS 保証会社利用初期費用カード決済通風良好

0 0 TV インターホン浴室乾燥機オートロック室内洗濯置システムキッチン住戸エレベーター洗面所独立コンロ宅配ボックス光ファイバー外壁タイル BS CS 防犯カメラ保証不要敷金保証金不要プール利用沿線利用徒歩 10分 24時間ゴミ敷地内ごみ置き場 5年都市ガス礼金

0 0 バルコニーエアコンクロゼットフローリング室内洗濯置シューズボックスシステムキッチン住戸外壁タイル入居ガスレンジ都市ガス

0 1 バストイレバルコニーエアコンクロゼットフローリング TV インターホン浴室乾燥機オートロック室内洗濯置シューズボックスシステムキッチン温水洗浄便座エレベーター洗面所独立コンロ駐輪場礼金不要バイク置場ガスレンジネット専用回線利用沿線利用徒歩 10分都市ガス BS

1 1 システム居室フローリングディンプルキーダブルロックキー 24時間換気システムフロントサービス耐火構造耐震構造利用沿線利用徒歩 5分 24時間ゴミ敷地内ごみ置き場セキュリティ会社加入済間接照明 LAN 保証会社利用初期費用カード決済通風良好

0 0 バストイレバルコニーフローリング TV インターホン浴室乾燥機オートロック室内洗濯置システムキッチン住戸エレベーター洗面所独立 CATV 入居防犯カメラグリル居室洋室仲エアコン納戸徒歩 10分専有面積南西

0 0 エアコンガスコンロ対応クロゼットフローリングシャワー洗面 TV インターホン浴室乾燥機室内洗濯置シューズボックス温水洗浄便座脱衣所洗面所独立駐輪場 CATV 保証不要 CATV インターネット平坦保証不要始発駅利用沿線利用徒歩 10分バス停徒歩 3分南西都市ガス BS

0 1 住戸浴室窓 24時間換気システムクロゼットルームシェア相談南面リビング学生相談一部フローリング利用沿線利用徒歩 5分駅徒歩 10分 24時間ゴミ敷地内ごみ置き場都市ガス洗面所ドア南面バルコニー BS 初期費用カード決済通風良好0 0 対面キッチン IH クッキングヒーターグリルウォークインクロゼット保証不要敷金二人入居相談ネット使用不要テラス保証不要 3年南面リビング浴室利用沿線利用バス停徒歩 3分敷地内ごみ置き場 LDK 都市ガス南面バルコニー天井シューズクロゼット礼金

0 0 バストイレバルコニーエアコンフローリング室内洗濯置シューズボックス焚機能浴室駐輪場宅配ボックス入居礼金不要敷金不要 IH クッキングヒーターエレベーター 2基利用沿線利用徒歩 5分都市ガス

1 1 コンロ宅配ボックス光ファイバー外壁タイル BS CS 防犯カメラ照明保証人不要敷金沿線利用ディンプルキー洗濯機耐火構造始発駅沿線利用徒歩 10分 24時間ゴミ敷地内ごみ置き場居室都市ガス高速ネット対応礼金保証会社利用初期費用カード決済

0 0 バルコニーエアコンクロゼットフローリングシューズボックス住戸エレベーター入居防犯カメラ分譲賃貸仲眺望良好保証不要利用沿線利用徒歩 5分敷地内ごみ置き場日勤管理敷金礼金不要通風良好

1 1 24時間換気システム南面リビング一部フローリング ISDN 対応耐火構造耐震構造利用徒歩 5分バス停徒歩 3分制震構造高層 24時間ゴミ敷地内ごみ置き場都市ガス南面バルコニー玄関収納 BS 高速ネット対応 LAN 初期費用カード決済

0 0 バイク置場南西住戸ガスレンジディンプルキー内装リフォーム済眺望良好照明センサー利用沿線利用徒歩 10分 24時間ゴミ敷地内ごみ置き場南西当社管理物件都市ガス洗面所ドア南面バルコニー室内物干礼金保証会社利用通風良好

0 0 駅平坦ネット使用不要ダブルロックキー保証不要 2年 24時間換気システム複層ガラス照明センサー利用徒歩 10分外壁サイディング 24時間ゴミ敷地内ごみ置き場南西都市ガス玄関収納年内入居年度内入居礼金初期費用カード決済

0 1 フローリング TV インターホン浴室乾燥機オートロック室内洗濯置陽良好シューズボックスシステムキッチン焚機能浴室住戸エレベーター洗面化粧台コンロ駐輪場宅配ボックス CATV 光ファイバー BS CS 対面キッチン沿線利用楽器相談ウッドデッキ利用徒歩 10分東南 LDK 居室都市ガス

LSTMの場合は、BLSTMに関する記述（1行）をネットワークから除外すれば良いので簡単に試すことができます。せっかくなので比較してみたところ、AUCは0.841でした。BLSTMは0.844だったので、若干上回る程度のようです。以前扱った記事ではデザイナーズマンション分類はAUCが0.83程度だったので、LSTMを使ったほうが精度が少しだけ高いようです。

追記

TensorBoardの検証のloglossを見てみます。

エポックに対してloglossが非常に荒れています。学習曲線としてはセオリー的にアウトなようです。一方で、検証のAUCはエポックに従って高まるようです。

AUCは良くなっているけどloglossが増え続けている。loglossが下がらないと過学習している可能性が高いので、これは過学習しているだけなのだろう。

参考文献

[R]ボージョレ・ヌーボーのコメントに対してLDATSパッケージを使って時系列トピックモデルを扱う

はじめに

先日、某勉強会でLTをしました。その際に10秒だけ紹介したRのパッケージについて記事を書いてみようと思います。

LDATSパッケージについて

時系列でのトピックモデルを推定することができるパッケージです。
やっていることとしてはLDAでトピックを推定して次元を減らし、そのトピックの多変量時系列に関してベイズ手法による変化点検知のためのパラメータ推定を行っているようです。GitHubの該当しそうなソースコードに多変量のデータに対するsoftmax関数での回帰をやっているとの記述がある。（multinomial Bayesian Time Series analysis）

元となっている論文を見る限り、BoW(Bag of Words)を想定して作っておらず、20~30程度のグループからなるデータに対して適用するのがちょうど良いです。アクセスログのページカテゴリや、マーケティングの顧客セグメントであればそんなに数は多くないので扱いやすいと思います。

データ

Webサイトから集めてきたボージョレ・ヌーボーのキャッチコピー14年分を今回は扱います。実は販売店側のキャッチコピーとワイン委員会が決めた評価が存在します。私の知っている世界は販売店側のキャッチコピーだけでした。

試してみた

今回はとにかく動くことだけを考えて、汚いコードとなっております。やっていることとしては、キャッチコピーを販売側とワイン委員会側のものを一つにつないで、数字を正規表現で「数字」に変換し、RMeCabで形態素解析をし、LDATS向けの形式のデータを作成していきます。
途中で、日本語の文字化け問題を回避するためにGoogle翻訳を使って単語名を置き換えています。
1時系列につき1文書となるようにデータを作っていく必要があるのですが、今回はボージョレ・ヌーボーのキャッチコピーなので最初から1時系列につき1文書となっているため都合が良いです。
データとソースコードはこちら。

library(tidyverse)
library(RMeCab)
library(LDATS)

wine_rating <- read_csv(file = "TimeseriesLDA/dataset.csv")

wine_rating <- wine_rating %>% filter(!is.na(sales_catch_copy_text),
                                      !is.na(ratings),
                                      year > 1999,
                                      !is.na(commission_text))

wine_rating$sales_catch_copy_text <- gsub(pattern = "[0-9]",
                                          replacement = "数字",
                                          x = wine_rating$sales_catch_copy_text)
wine_rating$commission_text <- gsub(pattern = "[0-9]",
                                    replacement = "数字",
                                    x = wine_rating$commission_text)

wine_rating$bind_text <- paste0(wine_rating$sales_catch_copy_text ,
                                wine_rating$commission_text) 


# Bag of wordsの生成
res <- docMatrixDF(wine_rating$bind_text,minFreq=3)
res <- data.frame(res)
# View(rownames(res))

# Google Spread Sheetの=GOOGLETRANSLATE(C18,"ja","en")で変換した英語のデータを読み込む
translate_df <- read_csv(file = "TimeseriesLDA/translate.csv",col_names = FALSE)
colnames(translate_df) <- c("word_ja", "word_en")

word_translate <- data.frame(word_ja=rownames(res))
word_translate <- word_translate %>% left_join(translate_df, by="word_ja")

rownames(res) <- word_translate$word_en

word_vector <- row.names(res)
colnames(res) <- wine_rating$year

# LDATSで扱えるデータ構造を作成
for (i in 1:nrow(res)) {
  nam <- paste( word_vector[i], sep = "")
  assign(nam, as.integer(res[i, ]))
}

# ここでの変数名が可視化の際に表示される
document_term_table <- data.frame(list(sa=sa,
                                       fruits=fruits,
                                       Greatness=Greatness,
                                       Can=Can,
                                       workmanship=workmanship,
                                       taste=taste,
                                       quality=quality,
                                       Year=Year,
                                       Thenumbers=Thenumbers,
                                       Highest=Highest,
                                       fruit=fruit,
                                       Great=Great,
                                       delicate=delicate,
                                       complexity=complexity,
                                       rich=rich,
                                       past=past,
                                       fragrance=fragrance
                                       ))

# 共変量データセット
document_covariate_table <- data.frame(list(year=as.integer(colnames(res))),
                                       list(rating=as.integer(wine_rating$ratings)))

test_set <- list(document_term_table=document_term_table,
                 document_covariate_table=document_covariate_table)

# 時系列トピックモデルの実行
r_LDATS <- LDA_TS(test_set,
                  topics = 3:6, 
                  nseeds = 2,
                  formulas = ~1,  
                  nchangepoints = 1:2,
                  timename = "year")
# 対数尤度などの出力
print(r_LDATS)

# 時系列トピックモデルの可視化（先行研究に準拠）
plot(r_LDATS)

library(tidyverse)

library(RMeCab)

library(LDATS)

wine_rating <- read_csv(file = "TimeseriesLDA/dataset.csv")

wine_rating <- wine_rating %>% filter(!is.na(sales_catch_copy_text),

!is.na(ratings),

year > 1999,

!is.na(commission_text))

wine_rating$sales_catch_copy_text <- gsub(pattern = "[0-9]",

replacement = "数字",

x = wine_rating$sales_catch_copy_text)

wine_rating$commission_text <- gsub(pattern = "[0-9]",

replacement = "数字",

x = wine_rating$commission_text)

wine_rating$bind_text <- paste0(wine_rating$sales_catch_copy_text ,

wine_rating$commission_text)

# Bag of wordsの生成

res <- docMatrixDF(wine_rating$bind_text,minFreq=3)

res <- data.frame(res)

# View(rownames(res))

# Google Spread Sheetの=GOOGLETRANSLATE(C18,"ja","en")で変換した英語のデータを読み込む

translate_df <- read_csv(file = "TimeseriesLDA/translate.csv",col_names = FALSE)

colnames(translate_df) <- c("word_ja", "word_en")

word_translate <- data.frame(word_ja=rownames(res))

word_translate <- word_translate %>% left_join(translate_df, by="word_ja")

rownames(res) <- word_translate$word_en

word_vector <- row.names(res)

colnames(res) <- wine_rating$year

# LDATSで扱えるデータ構造を作成

for (i in 1:nrow(res)) {

nam <- paste( word_vector[i], sep = "")

assign(nam, as.integer(res[i, ]))

}

# ここでの変数名が可視化の際に表示される

document_term_table <- data.frame(list(sa=sa,

fruits=fruits,

Greatness=Greatness,

Can=Can,

workmanship=workmanship,

taste=taste,

quality=quality,

Year=Year,

Thenumbers=Thenumbers,

Highest=Highest,

fruit=fruit,

Great=Great,

delicate=delicate,

complexity=complexity,

rich=rich,

past=past,

fragrance=fragrance

))

# 共変量データセット

document_covariate_table <- data.frame(list(year=as.integer(colnames(res))),

list(rating=as.integer(wine_rating$ratings)))

test_set <- list(document_term_table=document_term_table,

document_covariate_table=document_covariate_table)

# 時系列トピックモデルの実行

r_LDATS <- LDA_TS(test_set,

topics = 3:6,

nseeds = 2,

formulas = ~1,

nchangepoints = 1:2,

timename = "year")

# 対数尤度などの出力

print(r_LDATS)

# 時系列トピックモデルの可視化（先行研究に準拠）

plot(r_LDATS)

こちらは論文の図と同じものだとドキュメントの説明にあったので、論文の説明を見る限り、表すものとしては以下のようです。

一番上の積み上げグラフはトピックごとの単語の割合を表しています。
二番目の折れ線グラフはLDAによって推定されたトピックの時系列推移です。
三番目のヒストグラムは二番目の時系列における変化点を集計したものです。
四番目の折れ線グラフはモデルが推定したトピック割合の変化点の前後での推移です。

今回の図では文字が潰れていて見にくいですが、

トピック1はボキャブラリーが比較的リッチなコメント（「フルーティー」「フレグランス」「複雑」）
トピック2は数字を用いたコメント（「何年に一度の！」みたいな）
トピック3はボキャブラリーが貧相なコメント（「すごい！」みたいな）

のようです。
二番目の折れ線グラフを見る限り、周期的に数字を用いたコメントが現れているように思われます。四番目の折れ線グラフの変化点を見る限り、近年は数字を用いたコメントが相対的に減ってきて、リッチなボキャブラリーになってきているようです。

おわりに

時系列トピックモデルをカジュアルに試せる面白そうなパッケージだなと思い、LDATSパッケージを触ってみましたが、そもそもBoWなどを想定して作られているパッケージではないので、単語数が多いような分析ではそもそも可視化ができず使いにくいだろうなと思いました。マーケティングなどでユーザーのセグメントの推移を分析したい場合などにちょうど良いのだろうと思われます。

参考情報

[1] Long‐term community change through multiple rapid transitions in a desert rodent community
[2] Latent Dirichlet Allocation coupled with Bayesian Time Series analyses
[3] Package ‘LDATS’

R advent calendar 2019 RSelenium、jpmesh、sfパッケージで東京23区の事故物件を分析してみよう！

はじめに

今回で3回目となるR advent Calendarですが、前回は「一発屋芸人の検索トレンドのデータ」を扱い、前々回は「ポケモンのデータ」を扱いました。今回は人の命に関わるようなデータを扱ってみたいと思い、某サイトから東京都の23区内における事故物件の住所と詳細を集めてきました。どのようなエリアが事故が起きやすいのかの分析を行います。(以下では、事故物件をAP(Accident Property)と呼びます。)

R advent calendar用に向けたデータ収集コードができたが、データがけっこう汚いのでアノテーション工数を覚悟したほうが良さそう。
今回は人の命に関わるデータですね。

— SKUE (@Mr_Sakaue) December 8, 2019

分析工程

・データの収集
・データの整形
・可視化
・分析

データの収集

APに関する情報を某サイトより集める必要があります。そこで必要なライブラリとしては、RSeleniumやtwitteRがあげられます。
twitteRが必要な理由は、APに関するサイトにAPの一覧ページがなく、公式アカウントがAPに関するページのリンクを投稿しているところにあります。ただ、私が以前使っていた頃とはTwitterAPIの仕様が変わり、3ヶ月よりも前の情報にアクセスできなくなっていました。そのため、今後のデータに関してはTwitterAPIでいいのですが、過去のものに関しては別アプローチが必要となります。
また、APに関するサイトはJavaScriptで地図が表示されているだけなので、RSeleniumを使って地図をクリックさせ、表示された情報をスクレイピングするなどの処理が必要となります。
当初の想定ではTwitterのデータ収集、リンク先の情報をRSeleniumでスクレイピングするだけの簡単な仕事だと思っていたのですが、過去のデータにアクセスできないので、地図上で一つ一つ見つけていくためにRSeleniumだけで頑張ることにしました。(私の過去のアドカレ史上、一番面倒なデータとなりました。)

誰もやらないと思うのですが、一応手順を記しておきます。

RSeleniumだけでMapからAPの情報を抽出するための手順
1.都内の住所一覧を収集
2.検索窓に住所を入力
3.検索結果一覧の上位5件をクリック
4.一度地図を引くことでAPを広い範囲で捉えれるようにする
5.APのマークの要素を取得し、1件ずつクリックし、表示されたAPの情報をデータフレームに格納する

こちらが取得できたデータです。

RSeleniumでAPの
・住所
・発生時期（フリーテキスト）
・AP詳細
・投稿日
を集めることができるので、その住所データに対して、Yahoo!のジオコードに関するAPIを利用します。(利用申請が必要なはずです。4年前くらいに申請していたのでそのまま使えました。)
Yahoo!のAPIを使えば、住所から緯度経度の情報を取得することができます。

APの緯度経度がわかれば、jpmeshパッケージを用いて1kmメッシュやら10kmメッシュやらのメッシュデータに変換することができます。
jpmeshを用いてメッシュデータに変換し、メッシュ単位でAPの発生件数を集計します。

データ収集用のソースコードは思いのほか長くなってしまったので、GitHubにあげておきました。
https://github.com/KamonohashiPerry/r_advent_calendar_2019

ここで再度、手順を整理しておきましょう。

Twitterに出てきたものだけを取得（直近3ヶ月〜）する場合、
run_tweet_collect.RでTweetを収集
↓
run_selenium.RでAPの情報をスクレイピング
↓
run_map_api.Rで住所から緯度経度の取得
↓
making_mesh_data_and_download_other_data.Rで緯度経度からメッシュデータへの変換、その他の人口データや地価データと接続をします。

直近3ヶ月以前のものを取得する場合、
run_selenium_from_map.Rで地図上から直接APの情報を取得する
↓
making_mesh_data_and_download_other_data.Rで緯度経度からメッシュデータへの変換、その他の人口データや地価データと接続をします。

データの整形

APの1kmメッシュデータを手に入れたら、kokudosuuchiパッケージを使って国土地理院の収集したデータをつなぎこみます。手順としては、以下のとおりです。

まずは推計人口というそのエリアの人口の予測値です。今回は2010年のものを抽出しました。こちらは1kmメッシュのデータなので、変換することなく使えて都合が良いです。

library(kokudosuuchi)

# 推計人口
estimated_population <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/m1000/m1000-17/m1000-17_GML.zip",
                                      cache_dir = "cached_zip")

estimated_population_df <- data.frame(one_k_mesh=estimated_population$Mesh3_POP_00$MESH_ID,
                                         city_code=estimated_population$Mesh3_POP_00$CITY_CODE,
                                         population=estimated_population$Mesh3_POP_00$POP2010)

estimated_population_df <- estimated_population_df %>% 
                                mutate(one_k_mesh=as.integer(one_k_mesh))

library(kokudosuuchi)

# 推計人口

estimated_population <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/m1000/m1000-17/m1000-17_GML.zip",

cache_dir = "cached_zip")

estimated_population_df <- data.frame(one_k_mesh=estimated_population$Mesh3_POP_00$MESH_ID,

city_code=estimated_population$Mesh3_POP_00$CITY_CODE,

population=estimated_population$Mesh3_POP_00$POP2010)

estimated_population_df <- estimated_population_df %>%

mutate(one_k_mesh=as.integer(one_k_mesh))

続いて、2015年の公示地価を抽出しました。

# 公示地価
published_land_price <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/L01/L01-15/L01-15_GML.zip",
                                    cache_dir = "cached_zip")

published_land_price_df <- data.frame(price=as.integer(published_land_price$`L01-15`$L01_006),
                                      address=as.character(published_land_price$`L01-15`$L01_019),
                                      area_code=as.integer(published_land_price$`L01-15`$L01_017))

for (i in 1:nrow(published_land_price_df)) {
  published_land_price_df$latitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][2])
  published_land_price_df$longitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][1])
}

# 東京23区の地価に絞り込む
published_land_price_df_tokyo <- published_land_price_df %>% filter(stringr::str_detect(address,"東京"),
                                                                    stringr::str_detect(address,"区"))
published_land_price_df_tokyo$address_fixed <- gsub(pattern = "　", replacement = "", published_land_price_df_tokyo$address)
published_land_price_df_tokyo$address_fixed <- stringi::stri_trans_nfkc(published_land_price_df_tokyo$address_fixed)

for (i in 1:nrow(published_land_price_df_tokyo)) {
  published_land_price_df_tokyo$address_fixed[i] <- strsplit(published_land_price_df_tokyo$address_fixed,split="[0-9]")[[i]][1]
}

# 公示地価

published_land_price <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/L01/L01-15/L01-15_GML.zip",

cache_dir = "cached_zip")

published_land_price_df <- data.frame(price=as.integer(published_land_price$`L01-15`$L01_006),

address=as.character(published_land_price$`L01-15`$L01_019),

area_code=as.integer(published_land_price$`L01-15`$L01_017))

for (i in 1:nrow(published_land_price_df)) {

published_land_price_df$latitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][2])

published_land_price_df$longitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][1])

}

# 東京23区の地価に絞り込む

published_land_price_df_tokyo <- published_land_price_df %>% filter(stringr::str_detect(address,"東京"),

stringr::str_detect(address,"区"))

published_land_price_df_tokyo$address_fixed <- gsub(pattern = "　", replacement = "", published_land_price_df_tokyo$address)

published_land_price_df_tokyo$address_fixed <- stringi::stri_trans_nfkc(published_land_price_df_tokyo$address_fixed)

for (i in 1:nrow(published_land_price_df_tokyo)) {

published_land_price_df_tokyo$address_fixed[i] <- strsplit(published_land_price_df_tokyo$address_fixed,split="[0-9]")[[i]][1]

}

こちらはメッシュデータではないので、緯度経度の情報から1kmメッシュのデータに変換する必要があります。後で行います。

単純集計・可視化

今回のデータセットのデータ数は3919件です。本当は7000件以上はあると思われますが、マップから取ってくるという勝手上、なかなか全てを取り切ることができませんでした。

まずは、1kmメッシュごとのAP発生件数のヒストグラムです。

g <- ggplot(data = estimated_population_df_tokyo %>% filter(accident_count>0),
            aes(accident_count)) + theme_set(theme_bw(base_size = 14,base_family="HiraKakuProN-W3"))
g <- g +geom_histogram(bins = 50)
g <- g + ggtitle("1kmメッシュごとのAP件数のヒストグラム")
g

g <- ggplot(data = estimated_population_df_tokyo %>% filter(accident_count>0),

aes(accident_count)) + theme_set(theme_bw(base_size = 14,base_family="HiraKakuProN-W3"))

g <- g +geom_histogram(bins = 50)

g <- g + ggtitle("1kmメッシュごとのAP件数のヒストグラム")

1kmメッシュにおける人口のヒストグラムです。

公示地価のヒストグラムです。

1kmメッシュにおける人口あたりのAP件数のヒストグラムです。

分析

ここでは、色々な軸でAPのデータに向き合ってみようと思います。

APの発生件数の集計

人口が多いところがAPの発生件数が多いところだと思われますが、とりあえず確認します。

> estimated_population_df_tokyo %>% group_by(市区町村) %>% summarize(sum_accident=sum(accident_count), population=sum(population), mean_accident_percapita=mean(accident_percapita_100)) %>% arrange(desc(sum_accident))
Simple feature collection with 49 features and 4 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 139.15 ymin: 35.5 xmax: 139.9125 ymax: 35.84167
epsg (SRID):    NA
proj4string:    NA
# A tibble: 49 x 5
   市区町村 sum_accident population mean_accident_percap…                                                                     poligon
   <chr>           <dbl>      <dbl>                 <dbl>                                                                  <GEOMETRY>
 1 世田谷区          343    952572.                0.0337 POLYGON ((139.675 35.60833, 139.6875 35.60833, 139.6875 35.6, 139.675 35.6…
 2 新宿区            302    398145.                0.0795 POLYGON ((139.7 35.69167, 139.7 35.68333, 139.6875 35.68333, 139.6875 35.6…
 3 大田区            291    803602.                0.03   MULTIPOLYGON (((139.775 35.55, 139.7875 35.55, 139.7875 35.55833, 139.775 …
 4 板橋区            283    471726.                0.0570 POLYGON ((139.7125 35.75, 139.7125 35.74167, 139.7 35.74167, 139.7 35.75, …
 5 練馬区            276    575022.                0.0484 POLYGON ((139.6875 35.74167, 139.6875 35.73333, 139.675 35.73333, 139.6625…
 6 杉並区            268    549477.                0.0477 POLYGON ((139.65 35.675, 139.65 35.66667, 139.6375 35.66667, 139.6375 35.6…
 7 中野区            214    332012.                0.0619 MULTIPOLYGON (((139.675 35.69167, 139.675 35.68333, 139.6625 35.68333, 139…
 8 豊島区            203    278974.                0.0723 POLYGON ((139.725 35.725, 139.725 35.71667, 139.7125 35.71667, 139.7125 35…
 9 港区              159    188168.                0.137  MULTIPOLYGON (((139.7875 35.63333, 139.7875 35.625, 139.775 35.625, 139.77…
10 品川区            159    325433.                0.0382 POLYGON ((139.7625 35.59167, 139.7625 35.58333, 139.75 35.58333, 139.75 35…
# … with 39 more rows

> estimated_population_df_tokyo %>% group_by(市区町村) %>% summarize(sum_accident=sum(accident_count), population=sum(population), mean_accident_percapita=mean(accident_percapita_100)) %>% arrange(desc(sum_accident))

Simple feature collection with 49 features and 4 fields

geometry type: GEOMETRY

dimension: XY

bbox: xmin: 139.15 ymin: 35.5 xmax: 139.9125 ymax: 35.84167

epsg (SRID): NA

proj4string: NA

# A tibble: 49 x 5

市区町村 sum_accident population mean_accident_percap… poligon

1 世田谷区 343 952572. 0.0337 POLYGON ((139.675 35.60833, 139.6875 35.60833, 139.6875 35.6, 139.675 35.6…

2 新宿区 302 398145. 0.0795 POLYGON ((139.7 35.69167, 139.7 35.68333, 139.6875 35.68333, 139.6875 35.6…

3 大田区 291 803602. 0.03 MULTIPOLYGON (((139.775 35.55, 139.7875 35.55, 139.7875 35.55833, 139.775 …

4 板橋区 283 471726. 0.0570 POLYGON ((139.7125 35.75, 139.7125 35.74167, 139.7 35.74167, 139.7 35.75, …

5 練馬区 276 575022. 0.0484 POLYGON ((139.6875 35.74167, 139.6875 35.73333, 139.675 35.73333, 139.6625…

6 杉並区 268 549477. 0.0477 POLYGON ((139.65 35.675, 139.65 35.66667, 139.6375 35.66667, 139.6375 35.6…

7 中野区 214 332012. 0.0619 MULTIPOLYGON (((139.675 35.69167, 139.675 35.68333, 139.6625 35.68333, 139…

8 豊島区 203 278974. 0.0723 POLYGON ((139.725 35.725, 139.725 35.71667, 139.7125 35.71667, 139.7125 35…

9 港区 159 188168. 0.137 MULTIPOLYGON (((139.7875 35.63333, 139.7875 35.625, 139.775 35.625, 139.77…

10 品川区 159 325433. 0.0382 POLYGON ((139.7625 35.59167, 139.7625 35.58333, 139.75 35.58333, 139.75 35…

# … with 39 more rows

世田谷区は最も人口が多いことから、AP発生件数では一番となっています。続いて、歌舞伎町などがある新宿が来ています。しかしながら、人口に占めるAP発生件数で言うと、港区がかなり高く出ているのがわかります。

人口あたりのAP件数

ここでは、メッシュデータをsfパッケージ用のオブジェクトに変換して、1kmにおける人口あたりのAP発生割合を可視化しています。

library(jpmesh)
library(sf)
library(tidyverse)
library(kokudosuuchi)
library(jpndistrict)
library(mapview)

# load dataset
load(file = "accident_df_with_coordinate_from_map.RData")

# メッシュデータの追加
accident_df <- accident_df %>% mutate(k_mesh="")
accident_df <- accident_df %>% filter(latitude != "", !is.na(upload_date))

# 緯度、軽度の情報から任意のメッシュのデータを得る
for (i in 1:nrow(accident_df)) {
  accident_df$k_mesh[i] <- jpmesh::coords_to_mesh(
                                        as.numeric(accident_df$longitude[i]),
                                        as.numeric(accident_df$latitude[i]),
                                        mesh_size = "1km")
}

# メッシュを整数にする
accident_df <- accident_df %>% mutate(k_mesh=as.integer(k_mesh))

# 東京都に絞る
accident_df <- accident_df %>% filter(stringr::str_detect(formatted_address,"東京"))
accident_df %>% group_by(k_mesh) %>% count() %>% arrange(desc(n))

# 公示地価
published_land_price <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/L01/L01-15/L01-15_GML.zip",
                                    cache_dir = "cached_zip")

published_land_price_df <- data.frame(price=as.integer(published_land_price$`L01-15`$L01_006),
                                      address=as.character(published_land_price$`L01-15`$L01_019),
                                      area_code=as.integer(published_land_price$`L01-15`$L01_017))

for (i in 1:nrow(published_land_price_df)) {
  published_land_price_df$latitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][2])
  published_land_price_df$longitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][1])
}

# 東京の地価に絞り込む
published_land_price_df_tokyo <- published_land_price_df %>% filter(stringr::str_detect(address,"東京"),
                                                                    stringr::str_detect(address,"区"))
published_land_price_df_tokyo$address_fixed <- gsub(pattern = "　", replacement = "", published_land_price_df_tokyo$address)
published_land_price_df_tokyo$address_fixed <- stringi::stri_trans_nfkc(published_land_price_df_tokyo$address_fixed)

for (i in 1:nrow(published_land_price_df_tokyo)) {
  published_land_price_df_tokyo$address_fixed[i] <- strsplit(published_land_price_df_tokyo$address_fixed,split="[0-9]")[[i]][1]
}

# 緯度経度を1kmメッシュにする
for (i in 1:nrow(published_land_price_df)) {
  published_land_price_df$mesh[i] <- jpmesh::coords_to_mesh(
    as.numeric(published_land_price_df$longitude[i]),
    as.numeric(published_land_price_df$latitude[i]),
    mesh_size = "1km")
}

published_land_price_df <- published_land_price_df %>% mutate(mesh=as.integer(mesh))
published_land_price_df_summary <- published_land_price_df %>% group_by(mesh) %>% summarise(mean_price=mean(price))

# 推計人口
estimated_population <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/m1000/m1000-17/m1000-17_GML.zip",
                                   cache_dir = "cached_zip")

estimated_population_df <- data.frame(one_k_mesh=estimated_population$Mesh3_POP_00$MESH_ID,
                                       city_code=estimated_population$Mesh3_POP_00$CITY_CODE,
                                       population=estimated_population$Mesh3_POP_00$POP2010)

estimated_population_df <- estimated_population_df %>% 
                                mutate(one_k_mesh=as.integer(one_k_mesh))

accident_summary <- accident_df %>%
                        group_by(k_mesh) %>%
                        summarise(accident_count=n()) %>% 
                        arrange(desc(accident_count))

estimated_population_df <- estimated_population_df %>% 
                                        left_join(accident_summary,
                                                  by = c("one_k_mesh"="k_mesh"))

estimated_population_df$accident_count <- replace_na(estimated_population_df$accident_count, 0)

estimated_population_df <- estimated_population_df %>% mutate(accident_percapita=accident_count/population)

estimated_population_df <- estimated_population_df %>% 
                              left_join(published_land_price_df_summary,
                                        by=c("one_k_mesh"="mesh"))

estimated_population_df <- estimated_population_df %>% 
                                mutate(accident_ratio_class=if_else(estimated_population_df$accident_percapita > 0.05, "0.05以上",
                                                                    if_else(estimated_population_df$accident_percapita > 0.01, "0.01以上","0.01未満")))
estimated_population_df <- estimated_population_df %>% mutate(accident_percapita_100 = round(estimated_population_df$accident_percapita*100,2)  )


estimated_population_df_tokyo <- estimated_population_df %>% filter(city_code>=13000, city_code<13300)
estimated_population_df_tokyo <- estimated_population_df_tokyo %>% mutate(poligon="")

for (i in 1:nrow(estimated_population_df_tokyo)) {
  estimated_population_df_tokyo$poligon[i] <- export_mesh(estimated_population_df_tokyo$one_k_mesh[i])
}
# sfオブジェクトに変換
estimated_population_df_tokyo <- sf::st_as_sf(estimated_population_df_tokyo)

# マップでの可視化
estimated_population_df_tokyo %>% mapview::mapview(zcol = "accident_percapita_100")

100

101

102

103

104

105

106

107

108

library(jpmesh)

library(sf)

library(tidyverse)

library(kokudosuuchi)

library(jpndistrict)

library(mapview)

# load dataset

load(file = "accident_df_with_coordinate_from_map.RData")

# メッシュデータの追加

accident_df <- accident_df %>% mutate(k_mesh="")

accident_df <- accident_df %>% filter(latitude != "", !is.na(upload_date))

# 緯度、軽度の情報から任意のメッシュのデータを得る

for (i in 1:nrow(accident_df)) {

accident_df$k_mesh[i] <- jpmesh::coords_to_mesh(

as.numeric(accident_df$longitude[i]),

as.numeric(accident_df$latitude[i]),

mesh_size = "1km")

}

# メッシュを整数にする

accident_df <- accident_df %>% mutate(k_mesh=as.integer(k_mesh))

# 東京都に絞る

accident_df <- accident_df %>% filter(stringr::str_detect(formatted_address,"東京"))

accident_df %>% group_by(k_mesh) %>% count() %>% arrange(desc(n))

# 公示地価

published_land_price <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/L01/L01-15/L01-15_GML.zip",

cache_dir = "cached_zip")

published_land_price_df <- data.frame(price=as.integer(published_land_price$`L01-15`$L01_006),

address=as.character(published_land_price$`L01-15`$L01_019),

area_code=as.integer(published_land_price$`L01-15`$L01_017))

for (i in 1:nrow(published_land_price_df)) {

published_land_price_df$latitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][2])

published_land_price_df$longitude[i] <- as.character(published_land_price$`L01-15`$geometry[[i]][1])

}

# 東京の地価に絞り込む

published_land_price_df_tokyo <- published_land_price_df %>% filter(stringr::str_detect(address,"東京"),

stringr::str_detect(address,"区"))

published_land_price_df_tokyo$address_fixed <- gsub(pattern = "　", replacement = "", published_land_price_df_tokyo$address)

published_land_price_df_tokyo$address_fixed <- stringi::stri_trans_nfkc(published_land_price_df_tokyo$address_fixed)

for (i in 1:nrow(published_land_price_df_tokyo)) {

published_land_price_df_tokyo$address_fixed[i] <- strsplit(published_land_price_df_tokyo$address_fixed,split="[0-9]")[[i]][1]

}

# 緯度経度を1kmメッシュにする

for (i in 1:nrow(published_land_price_df)) {

published_land_price_df$mesh[i] <- jpmesh::coords_to_mesh(

as.numeric(published_land_price_df$longitude[i]),

as.numeric(published_land_price_df$latitude[i]),

mesh_size = "1km")

}

published_land_price_df <- published_land_price_df %>% mutate(mesh=as.integer(mesh))

published_land_price_df_summary <- published_land_price_df %>% group_by(mesh) %>% summarise(mean_price=mean(price))

# 推計人口

estimated_population <- getKSJData("http://nlftp.mlit.go.jp/ksj/gml/data/m1000/m1000-17/m1000-17_GML.zip",

cache_dir = "cached_zip")

estimated_population_df <- data.frame(one_k_mesh=estimated_population$Mesh3_POP_00$MESH_ID,

city_code=estimated_population$Mesh3_POP_00$CITY_CODE,

population=estimated_population$Mesh3_POP_00$POP2010)

estimated_population_df <- estimated_population_df %>%

mutate(one_k_mesh=as.integer(one_k_mesh))

accident_summary <- accident_df %>%

group_by(k_mesh) %>%

summarise(accident_count=n()) %>%

arrange(desc(accident_count))

estimated_population_df <- estimated_population_df %>%

left_join(accident_summary,

by = c("one_k_mesh"="k_mesh"))

estimated_population_df$accident_count <- replace_na(estimated_population_df$accident_count, 0)

estimated_population_df <- estimated_population_df %>% mutate(accident_percapita=accident_count/population)

estimated_population_df <- estimated_population_df %>%

left_join(published_land_price_df_summary,

by=c("one_k_mesh"="mesh"))

estimated_population_df <- estimated_population_df %>%

mutate(accident_ratio_class=if_else(estimated_population_df$accident_percapita > 0.05, "0.05以上",

if_else(estimated_population_df$accident_percapita > 0.01, "0.01以上","0.01未満")))

estimated_population_df <- estimated_population_df %>% mutate(accident_percapita_100 = round(estimated_population_df$accident_percapita*100,2) )

estimated_population_df_tokyo <- estimated_population_df %>% filter(city_code>=13000, city_code<13300)

estimated_population_df_tokyo <- estimated_population_df_tokyo %>% mutate(poligon="")

for (i in 1:nrow(estimated_population_df_tokyo)) {

estimated_population_df_tokyo$poligon[i] <- export_mesh(estimated_population_df_tokyo$one_k_mesh[i])

}

# sfオブジェクトに変換

estimated_population_df_tokyo <- sf::st_as_sf(estimated_population_df_tokyo)

# マップでの可視化

estimated_population_df_tokyo %>% mapview::mapview(zcol = "accident_percapita_100")

こちらのmapviewパッケージで作ったマップはインタラクティブにいじることができます。ぜひ関心のあるエリアでいじってみてください。

1kmメッシュ人口あたりのAP発生件数（×100）の可視化

<br />

1kmメッシュでのAP発生件数の可視化

<br />

比率ベースで、色の明るいメッシュのところを見ると、港区、中央区、新宿区、渋谷区などがAPが発生しやすいようです。件数ベースで言うと新宿が一番多いですね。
一番色が明るい港区はてっきり六本木ではないかと思ったのですが、新橋から日比谷にかけたエリアでした。会社員による自○が多いようです。恐ろしいものです。

APの名前の集計

APの名前を集計してみます。これは別にこの名前だからAPになりやすいというわけではなく、単純に数が多いだけの可能性がありますし、実際にそうだろうと思われます。AP発生率を知るには、APではないものも含めた全物件名に占めるAP発生件名を手に入れないといけませんが、全物件名を収集するのが難しいことから単純に頻度の集計となります。今回は、wordcloud2パッケージを使って、ワードクラウドにしてみます。文字が大きいと頻度が高いものとなります。

library('wordcloud2')
library(RMeCab)

load(file = "accident_df_with_coordinate__from_map.RData")

address_bow <- docMatrixDF(accident_df$address,minFreq=1)
word_list <- data.frame(word=rownames(address_bow), count=rowSums(address_bow))
word_list %>% filter(count<100, !stringr::str_detect(word, "[0-9]")) %>% wordcloud2()

library('wordcloud2')

library(RMeCab)

load(file = "accident_df_with_coordinate__from_map.RData")

address_bow <- docMatrixDF(accident_df$address,minFreq=1)

word_list <- data.frame(word=rownames(address_bow), count=rowSums(address_bow))

word_list %>% filter(count<100, !stringr::str_detect(word, "[0-9]")) %>% wordcloud2()

ハイツ、荘、コーポ、マンション、号棟、アパート、ハウスなどが多く出現しているようです。ただ、物件の名前としても頻度が高いとも考えられますね。

地価と人口あたりのAP発生件数の関係

ここでは地価のデータとAP発生の関係性について見てみます。

estimated_population_df_tokyo <- estimated_population_df_tokyo %>% mutate(bin = ntile(mean_price, 10))
vis_df <- estimated_population_df_tokyo %>% filter(accident_count>0) %>% group_by(bin) %>% summarize(mean_accident_percapita=mean(accident_percapita),
                                                                                                     count=n(),
                                                                                                     mean_price=mean(mean_price))

g <- ggplot(data = vis_df, aes(x =  bin, y =mean_accident_percapita)) + theme_set(theme_bw(base_size = 14,base_family="HiraKakuProN-W3"))
g <- g +geom_point() + stat_smooth(se=T,fullrange = T,level = 0.95)
g <- g + ggtitle("地価の階級値と平均人口あたりのAP発生件数")
g

estimated_population_df_tokyo <- estimated_population_df_tokyo %>% mutate(bin = ntile(mean_price, 10))

vis_df <- estimated_population_df_tokyo %>% filter(accident_count>0) %>% group_by(bin) %>% summarize(mean_accident_percapita=mean(accident_percapita),

count=n(),

mean_price=mean(mean_price))

g <- ggplot(data = vis_df, aes(x = bin, y =mean_accident_percapita)) + theme_set(theme_bw(base_size = 14,base_family="HiraKakuProN-W3"))

g <- g +geom_point() + stat_smooth(se=T,fullrange = T,level = 0.95)

g <- g + ggtitle("地価の階級値と平均人口あたりのAP発生件数")

地価の階級値（10個のパーセンタイルに分割）を横軸に、縦軸に人口あたりAP発生数をおくと、地価が上がるに従い人口あたりAP発生数が高まる傾向があります。これは、人口密度が高く地価の高いところではAPが発生しやすいということを示しているのではないでしょうか。人口密度が高いと地価があがる、人口密度が高いと治安が悪くなるという可能性が考えられます。

APの詳細の集計

ここではAPになってしまった詳細の内容について先ほどと同様に形態素解析を行いワードクラウドにしてみます。

accident_bow <- docMatrixDF(accident_df$detail,minFreq=1)
word_list <- data.frame(word=rownames(accident_bow), count=rowSums(accident_bow))
word_list %>% filter(count<300, !stringr::str_detect(word, "[0-9]")) %>% wordcloud2()

accident_bow <- docMatrixDF(accident_df$detail,minFreq=1)

word_list <- data.frame(word=rownames(accident_bow), count=rowSums(accident_bow))

word_list %>% filter(count<300, !stringr::str_detect(word, "[0-9]")) %>% wordcloud2()

どうやら孤独死が多いようです。高齢者の人口構成比が関係しているのだろうと思われます。

APの発生時期に関するテキストマイニング

ここでは、発生時期に含まれる四桁の数字を集計して、何年くらいのAPが多いのかをざっくりと掴みます。

date_bow <- docMatrixDF(accident_df$date,minFreq=1)
word_list <- data.frame(word=rownames(date_bow), count=rowSums(date_bow))
word_list %>% filter(stringr::str_detect(word, "[0-9][0-9][0-9][0-9]")) %>% wordcloud2()

date_bow <- docMatrixDF(accident_df$date,minFreq=1)

word_list <- data.frame(word=rownames(date_bow), count=rowSums(date_bow))

word_list %>% filter(stringr::str_detect(word, "[0-9][0-9][0-9][0-9]")) %>% wordcloud2()

どうやら昔のデータはあまり登録されていないようです。記憶が確かではないかもしれませんし、古すぎるものは消されているのかもしれませんね。あのサイトはユーザー生成コンテンツ（UGC）なので、投稿する人はそこまで昔のことをわざわざ投稿しないのかもしれないですね。

APの詳細に関するテキストマイニング

ここではトピック数10として、topicmodelsパッケージを使いLDAを行います。

library(RMeCab)
library(topicmodels)

accident_df$detail_fixed <- gsub(x = accident_df$detail,
                                 pattern = "\\(|\\)|（|）|/|\\.|:|_|-|[0-9]|[０-９]",
                                 replacement = " ")

# Bag of wordsの生成
res <- docMatrixDF(accident_df$detail_fixed,minFreq=4)
colnames(res) <- accident_df$k_mesh
res <- t(res)
res <- res[rowSums(res) >=1, ]

# トピック数
k <- 10
# 出力単語数
i <- 5
# 出力トピック数
j <- 5

# LDAの推定
LDA_estimate <- LDA(res, k,method="Gibbs",control=list(verbose=1))

# トピックごとの単語の出力
terms_each_topics <- data.frame(terms(LDA_estimate,i))

library(RMeCab)

library(topicmodels)

accident_df$detail_fixed <- gsub(x = accident_df$detail,

pattern = "\$|\$|（|）|/|\\.|:|_|-|[0-9]|[０-９]",

replacement = " ")

# Bag of wordsの生成

res <- docMatrixDF(accident_df$detail_fixed,minFreq=4)

colnames(res) <- accident_df$k_mesh

res <- t(res)

res <- res[rowSums(res) >=1, ]

# トピック数

k <- 10

# 出力単語数

i <- 5

# 出力トピック数

j <- 5

# LDAの推定

LDA_estimate <- LDA(res, k,method="Gibbs",control=list(verbose=1))

# トピックごとの単語の出力

terms_each_topics <- data.frame(terms(LDA_estimate,i))

なかなか恐ろしいキーワードが多いですが、なんとなくですがうまく分類されているのではないかと思われます。

トピック1は男性の不幸
トピック2は不動産屋に言われた告知事項
トピック3は孤独死
トピック4は病死
トピック5は火災・転落・事故
トピック6は事故のあった建物に関する記載
トピック7は腐乱した事例
トピック8は建物に関して不明であることの記載
トピック9は心理的瑕疵あり
トピック10は自○

となっているように思われます。まさかこのようなデータにトピックモデルを使うことになるとは。

おわりに

今回はR言語のみを用いて、APに関するデータを収集し、地図にプロットしたり他のメッシュデータとつなぎ合わせて分析をするなどしました。APが発生しやすいエリア、APと地価との関係、APのテキストマイニングなど興味深い結果が得られたと思います。
一つ残念なのは、時系列情報がフリーテキストなので、APがどのエリアでどの頻度で発生していくのかの分析のコストが高く、今回は時系列情報を用いた分析にチャレンジできませんでした。
今後はタクシーの需要推定の分析で行われているように、メッシュ単位でのAP発生確率の推定などを機械学習で行えると面白いなと思います。どなたか一緒にアノテーションしましょう！

それでは、どうか良い年末をお過ごし下さい！
メリークリスマス！

参考情報

数多くの方々の記事を見てどうにか仕上げることができました。感謝します。

[1]【追記あり】sfパッケージでシェープファイルを読み込んでmapviewパッケージで可視化するまで
[2]How to use mesh cord in R
[3]Rを使ってワードクラウドを作ってみました
[4]国土数値情報ダウンロードサービスWeb APIからデータを取得するためのRパッケージです
[5]東京の地価公示データを眺める
[6]Chapter 1 Introduction to spatial data in R
[7][翻訳] RSelenium vignette: RSeleniumの基本
[8]RからYahoo!のジオコーディングを利用する方法
[9]EMBEDDING A LEAFLET MAP ON WORDPRESS
[10]mapview advanced controls
[11]RSeleniumでChromeからファイルをダウンロードするディレクトリを指定する方法
[12]Selenium Serverが立ち上がらないときはportが被っているかも！？
[13]brew install selenium-server-standalone
[14]ナウでヤングなRの環境変数管理方法
[15]タクシードライバー向け需要予測について
[16]LDA with topicmodels package for R, how do I get the topic probability for each term?
[17]dplyr — 高速data.frame処理

[Python]機械学習などでテキストデータを特徴量にする際のソースコード集

テキストデータの特徴量化について

仕事ではテキストデータを多用するので、機械学習などで扱うためにテキストデータを特徴量にするためのアプローチを色々と整理してソースコードを残しておきたいと思います。今回はあくまでも私の知っているものだけなので、網羅性はないかもしれませんが悪しからず。
（2019/08/18 追記）Stackingをカジュアルに行えるvecstackというモジュールを用いた予測も試してみました。下の方の追記をご覧ください。

アプローチ

テキストデータを特徴量にする際のアプローチとしては、以下の3つが良く使っているものとなります。
・単語ベース
・クラスタ、トピック、分散表現ベース
・文書間の類似度ベース

今回扱うデータ

ひょんなことから、昨年10月くらいに取りためたマンションの施設情報のテキストです。

緑色が印象的な某不動産紹介サイトをクローリングしました。全部で1864件ほどの文書数となります。

加えて、デザイナーズマンションかどうかのフラグを作成しました（17%くらいがデザイナーズマンションの割合）。これでもって、マンションの施設情報からデザイナーズマンションかどうかを分類できるかチャレンジしたいと思います。
ここにデータを置いていますので、興味のある方はご利用ください。

今回扱うモデル

ランダムフォレストです。10foldsクロスバリデーションによるAUCの結果を各手法のスコアとして扱います。

こちらは、任意の手法に関して10foldsクロスバリデーションで実行し、AUCのグラフを生成してくれるソースコードです。主にscikit-learnのサイトに載っているものです。引数のclassifierをsklearnの任意のモデルのインスタンスで渡せば動きます。

def classification_and_roc_analysis(k, classifier, X, y):
    import numpy as np
    from scipy import interp
    import matplotlib.pyplot as plt
    
    from pylab import rcParams
    rcParams['figure.figsize'] = 8,8

    from sklearn.model_selection import StratifiedKFold
    from sklearn.metrics import roc_curve, auc
    from sklearn.ensemble import RandomForestClassifier
    
    # Classification and ROC analysis
    # Run classifier with cross-validation and plot ROC curves
    cv = StratifiedKFold(n_splits=k, random_state = 123)
    classifier = classifier

    tprs = []
    aucs = []
    mean_fpr = np.linspace(0, 1, 100)

    i = 0
    for train, test in cv.split(X, y):
        probas_ = classifier.fit(X[train], y[train]).predict_proba(X[test])
        # Compute ROC curve and area the curve
        fpr, tpr, thresholds = roc_curve(y[test], probas_[:, 1])
        tprs.append(interp(mean_fpr, fpr, tpr))
        tprs[-1][0] = 0.0
        roc_auc = auc(fpr, tpr)
        aucs.append(roc_auc)
        plt.plot(fpr, tpr, lw=1, alpha=0.3,
                 label='ROC fold %d (AUC = %0.2f)' % (i, roc_auc))

        i += 1
    plt.plot([0, 1], [0, 1], linestyle='--', lw=2, color='r',
             label='Chance', alpha=.8)

    mean_tpr = np.mean(tprs, axis=0)
    mean_tpr[-1] = 1.0
    mean_auc = auc(mean_fpr, mean_tpr)
    std_auc = np.std(aucs)
    plt.plot(mean_fpr, mean_tpr, color='b',
             label=r'Mean ROC (AUC = %0.2f $\pm$ %0.2f)' % (mean_auc, std_auc),
             lw=2, alpha=.8)

    std_tpr = np.std(tprs, axis=0)
    tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
    tprs_lower = np.maximum(mean_tpr - std_tpr, 0)
    plt.fill_between(mean_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2,
                     label=r'$\pm$ 1 std. dev.')

    plt.xlim([-0.05, 1.05])
    plt.ylim([-0.05, 1.05])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('Receiver operating characteristic example')
    plt.legend(loc="lower right")
    plt.show()

def classification_and_roc_analysis(k, classifier, X, y):

import numpy as np

from scipy import interp

import matplotlib.pyplot as plt

from pylab import rcParams

rcParams['figure.figsize'] = 8,8

from sklearn.model_selection import StratifiedKFold

from sklearn.metrics import roc_curve, auc

from sklearn.ensemble import RandomForestClassifier

# Classification and ROC analysis

# Run classifier with cross-validation and plot ROC curves

cv = StratifiedKFold(n_splits=k, random_state = 123)

classifier = classifier

tprs = []

aucs = []

mean_fpr = np.linspace(0, 1, 100)

i = 0

for train, test in cv.split(X, y):

probas_ = classifier.fit(X[train], y[train]).predict_proba(X[test])

# Compute ROC curve and area the curve

fpr, tpr, thresholds = roc_curve(y[test], probas_[:, 1])

tprs.append(interp(mean_fpr, fpr, tpr))

tprs[-1][0] = 0.0

roc_auc = auc(fpr, tpr)

aucs.append(roc_auc)

plt.plot(fpr, tpr, lw=1, alpha=0.3,

label='ROC fold %d (AUC = %0.2f)' % (i, roc_auc))

i += 1

plt.plot([0, 1], [0, 1], linestyle='--', lw=2, color='r',

label='Chance', alpha=.8)

mean_tpr = np.mean(tprs, axis=0)

mean_tpr[-1] = 1.0

mean_auc = auc(mean_fpr, mean_tpr)

std_auc = np.std(aucs)

plt.plot(mean_fpr, mean_tpr, color='b',

label=r'Mean ROC (AUC = %0.2f $\pm$ %0.2f)' % (mean_auc, std_auc),

lw=2, alpha=.8)

std_tpr = np.std(tprs, axis=0)

tprs_upper = np.minimum(mean_tpr + std_tpr, 1)

tprs_lower = np.maximum(mean_tpr - std_tpr, 0)

plt.fill_between(mean_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2,

label=r'$\pm$ 1 std. dev.')

plt.xlim([-0.05, 1.05])

plt.ylim([-0.05, 1.05])

plt.xlabel('False Positive Rate')

plt.ylabel('True Positive Rate')

plt.title('Receiver operating characteristic example')

plt.legend(loc="lower right")

plt.show()

単語ベース

シンプルに単語をそのまま特徴量にするというものですが、文書によっては単語数が多すぎて収集がつかないと思います。そこで単語を簡単に選択できるDocumentFeatureSelectionというパッケージを利用します。

# 事前にインストールしておく
pip install DocumentFeatureSelection
pip install JapaneseTokenizer
pip install neologdn

# 事前にインストールしておく

pip install DocumentFeatureSelection

pip install JapaneseTokenizer

pip install neologdn

このパッケージでは
・TF-IDFベースの特徴量選択
・PMI（Pointwise Mutual Information）ベースの特徴量選択
・SOA（Strength of association）ベースの特徴量選択
・BNS（Bi-Normal Separation）ベースの特徴量選択
を行うことができます。

まずは今回のベースラインとして、単語のカウントベースでの特徴量を扱いたいと思います。
その前に、GitHubに上がっているデータに対して以下のように簡単な前処理をしておきます。

%matplotlib inline
import pandas as pd
import numpy as np
import seaborn as sns
import copy
import warnings
warnings.filterwarnings('ignore')

# Data Import
corpus_data = pd.read_csv("designers_apartment.csv")
corpus_data = corpus_data[~corpus_data.text.isna()].reset_index(drop=True)

# 形態素解析
import collections
import MeCab
import mojimoji
from string import digits

remove_digits = str.maketrans('', '', digits)

tagger = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/")

def nouns_extract(line):
    keyword=[]
    node = tagger.parseToNode(line).next
    while node:
        if node.feature.split(",")[0] == "名詞":
            keyword.append(node.surface)
        node = node.next
    keyword = str(keyword).replace("', '"," ")
    keyword = keyword.replace("\'","")
    keyword = keyword.replace("[","")
    keyword = keyword.replace("]","")
    return keyword

# 全角を半角にする
corpus_data["text_fixed"] =  list(map(lambda text: mojimoji.zen_to_han(text, kana=False) , corpus_data.text))

# 数字を除外する
corpus_data["text_fixed"] =  list(map(lambda text: text.translate(remove_digits) , corpus_data.text_fixed))

# デザイナーズという表現を除外する
corpus_data["text_fixed"] =  corpus_data.text_fixed.str.replace("デザイナーズ", "")

# 形態素解析する
corpus_data["text_wakati"] = list(map(lambda text:nouns_extract(text) , list(corpus_data.text_fixed)))

%matplotlib inline

import pandas as pd

import numpy as np

import seaborn as sns

import copy

import warnings

warnings.filterwarnings('ignore')

# Data Import

corpus_data = pd.read_csv("designers_apartment.csv")

corpus_data = corpus_data[~corpus_data.text.isna()].reset_index(drop=True)

# 形態素解析

import collections

import MeCab

import mojimoji

from string import digits

remove_digits = str.maketrans('', '', digits)

tagger = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/")

def nouns_extract(line):

keyword=[]

node = tagger.parseToNode(line).next

while node:

if node.feature.split(",")[0] == "名詞":

keyword.append(node.surface)

node = node.next

keyword = str(keyword).replace("', '"," ")

keyword = keyword.replace("\'","")

keyword = keyword.replace("[","")

keyword = keyword.replace("]","")

return keyword

# 全角を半角にする

corpus_data["text_fixed"] = list(map(lambda text: mojimoji.zen_to_han(text, kana=False) , corpus_data.text))

# 数字を除外する

corpus_data["text_fixed"] = list(map(lambda text: text.translate(remove_digits) , corpus_data.text_fixed))

# デザイナーズという表現を除外する

corpus_data["text_fixed"] = corpus_data.text_fixed.str.replace("デザイナーズ", "")

# 形態素解析する

corpus_data["text_wakati"] = list(map(lambda text:nouns_extract(text) , list(corpus_data.text_fixed)))

ようやくベースラインの予測となります。以下のコードを実行すると、ROCが描かれた図がJupyter上で表示されます。

from sklearn.feature_extraction.text import CountVectorizer

# word count feature
countvector = CountVectorizer(min_df=0.01, max_df=0.90)
matrix = countvector.fit_transform(corpus_data.text_wakati)
countvector.get_feature_names()[:10]

X = matrix
y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

from sklearn.feature_extraction.text import CountVectorizer

# word count feature

countvector = CountVectorizer(min_df=0.01, max_df=0.90)

matrix = countvector.fit_transform(corpus_data.text_wakati)

countvector.get_feature_names()[:10]

X = matrix

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

AUC82%というのはベースラインとしてはなかなか強敵なのではないでしょうか。

さて、本題の特徴量選択パッケージの適用をするためのソースコードを以下に記します。

from scipy.sparse import csr_matrix
import JapaneseTokenizer
from DocumentFeatureSelection import interface

# `mecab-config` コマンドが通っているパスを書きます。場所がわからない時は`which mecab-config`で探してみてください。
path_mecab_config='/usr/local/bin'

# 辞書タイプを選びます。"neologd", "all", "ipaddic", "user", ""が選べます。
dictType = "neologd"

pos_condition = [('名詞',)]

mecab_wrapper = JapaneseTokenizer.MecabWrapper(
    dictType=dictType, 
    path_mecab_config=path_mecab_config
)

def tokenize_and_filtering(docs):
    assert isinstance(docs, list)
    tokenized_docs = [
        mecab_wrapper.tokenize(
            sentence=doc.strip('\n'),
            return_list=False
            )
        for doc in docs
    ]
    
    filtered_tokens = [
        mecab_wrapper.filter(
            parsed_sentence=tokenized_obj,
            pos_condition=pos_condition
        ).convert_list_object()
        for tokenized_obj in tokenized_docs
    ]

    return filtered_tokens

# デザイナーズマンションとそうでないマンションでテキストを分ける
processed_designer = tokenize_and_filtering(corpus_data.query(' designer_flag ==1 ').text_fixed.tolist())
processed_normal = tokenize_and_filtering(corpus_data.query(' designer_flag ==0 ').text_fixed.tolist())

input_labeled_docs_dict = {
    'designer': processed_designer,
    'normal': processed_normal
}

# TF-IDF base
tf_idf_scored_object = interface.run_feature_selection(
    input_dict=input_labeled_docs_dict,
    method='tf_idf',
    n_jobs=5
)
tf_idf_scored_df = pd.DataFrame(tf_idf_scored_object.ScoreMatrix2ScoreDictionary())

# PMI base
pmi_scored_object_cython = interface.run_feature_selection(
    input_dict=input_labeled_docs_dict,
    method='pmi',
    use_cython=True
)
pmi_scored_df = pd.DataFrame(pmi_scored_object_cython.ScoreMatrix2ScoreDictionary())

# SOA base
soa_scored_object_cython = interface.run_feature_selection(
    input_dict=input_labeled_docs_dict,
    method='soa',
    use_cython=True
)
soa_scored_df = pd.DataFrame(soa_scored_object_cython.ScoreMatrix2ScoreDictionary())

# BNS base
bns_scored_object = interface.run_feature_selection(
    input_dict=input_labeled_docs_dict,
    method='bns',
    n_jobs=5
)
bns_scored_df = pd.DataFrame(bns_scored_object.ScoreMatrix2ScoreDictionary())

from scipy.sparse import csr_matrix

import JapaneseTokenizer

from DocumentFeatureSelection import interface

# `mecab-config` コマンドが通っているパスを書きます。場所がわからない時は`which mecab-config`で探してみてください。

path_mecab_config='/usr/local/bin'

# 辞書タイプを選びます。"neologd", "all", "ipaddic", "user", ""が選べます。

dictType = "neologd"

pos_condition = [('名詞',)]

mecab_wrapper = JapaneseTokenizer.MecabWrapper(

dictType=dictType,

path_mecab_config=path_mecab_config

)

def tokenize_and_filtering(docs):

assert isinstance(docs, list)

tokenized_docs = [

mecab_wrapper.tokenize(

sentence=doc.strip('\n'),

return_list=False

)

for doc in docs

]

filtered_tokens = [

mecab_wrapper.filter(

parsed_sentence=tokenized_obj,

pos_condition=pos_condition

).convert_list_object()

for tokenized_obj in tokenized_docs

]

return filtered_tokens

# デザイナーズマンションとそうでないマンションでテキストを分ける

processed_designer = tokenize_and_filtering(corpus_data.query(' designer_flag ==1 ').text_fixed.tolist())

processed_normal = tokenize_and_filtering(corpus_data.query(' designer_flag ==0 ').text_fixed.tolist())

input_labeled_docs_dict = {

'designer': processed_designer,

'normal': processed_normal

}

# TF-IDF base

tf_idf_scored_object = interface.run_feature_selection(

input_dict=input_labeled_docs_dict,

method='tf_idf',

n_jobs=5

)

tf_idf_scored_df = pd.DataFrame(tf_idf_scored_object.ScoreMatrix2ScoreDictionary())

# PMI base

pmi_scored_object_cython = interface.run_feature_selection(

input_dict=input_labeled_docs_dict,

method='pmi',

use_cython=True

)

pmi_scored_df = pd.DataFrame(pmi_scored_object_cython.ScoreMatrix2ScoreDictionary())

# SOA base

soa_scored_object_cython = interface.run_feature_selection(

input_dict=input_labeled_docs_dict,

method='soa',

use_cython=True

)

soa_scored_df = pd.DataFrame(soa_scored_object_cython.ScoreMatrix2ScoreDictionary())

# BNS base

bns_scored_object = interface.run_feature_selection(

input_dict=input_labeled_docs_dict,

method='bns',

n_jobs=5

)

bns_scored_df = pd.DataFrame(bns_scored_object.ScoreMatrix2ScoreDictionary())

以上のソースコードを実行すれば、tf_idf_scored_df、pmi_scored_df、soa_scored_df、bns_scored_dfにスコアを付与された単語のリストが手に入ります。

ここでは各スコアに関してアドホックに閾値を設けて、特徴量として利用することにします。

TF-IDFベースの特徴量選択

tfidf_list = tf_idf_scored_df.query(' score > 0.001 ').feature.unique().tolist()
corpus_data["tfidf_base_list"] = list(map(lambda text: " ".join(list(set(text.split(" ") ) & set(tfidf_list))) , corpus_data.text_wakati))

# word count feature
countvector = CountVectorizer()
matrix = countvector.fit_transform(corpus_data.tfidf_base_list)

X = matrix
y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

tfidf_list = tf_idf_scored_df.query(' score > 0.001 ').feature.unique().tolist()

corpus_data["tfidf_base_list"] = list(map(lambda text: " ".join(list(set(text.split(" ") ) & set(tfidf_list))) , corpus_data.text_wakati))

# word count feature

countvector = CountVectorizer()

matrix = countvector.fit_transform(corpus_data.tfidf_base_list)

X = matrix

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

PMIベースの特徴量選択

pmi_list = pmi_scored_df.query(' score > 0.001 ').feature.unique().tolist()
corpus_data["pmi_base_list"] = list(map(lambda text: " ".join(list(set(text.split(" ") ) & set(pmi_list))) , corpus_data.text_wakati))

# word count feature
countvector = CountVectorizer()
matrix = countvector.fit_transform(corpus_data.pmi_base_list)

X = matrix
y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

pmi_list = pmi_scored_df.query(' score > 0.001 ').feature.unique().tolist()

corpus_data["pmi_base_list"] = list(map(lambda text: " ".join(list(set(text.split(" ") ) & set(pmi_list))) , corpus_data.text_wakati))

# word count feature

countvector = CountVectorizer()

matrix = countvector.fit_transform(corpus_data.pmi_base_list)

X = matrix

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

SOAベースの特徴量選択

soa_list = soa_scored_df.query(' score > 0.2 ').feature.unique().tolist()
corpus_data["soa_base_list"] = list(map(lambda text: " ".join(list(set(text.split(" ") ) & set(soa_list))) , corpus_data.text_wakati))

# word count feature
countvector = CountVectorizer()
matrix = countvector.fit_transform(corpus_data.soa_base_list)

X = matrix
y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

soa_list = soa_scored_df.query(' score > 0.2 ').feature.unique().tolist()

corpus_data["soa_base_list"] = list(map(lambda text: " ".join(list(set(text.split(" ") ) & set(soa_list))) , corpus_data.text_wakati))

# word count feature

countvector = CountVectorizer()

matrix = countvector.fit_transform(corpus_data.soa_base_list)

X = matrix

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

BNSベースの特徴量選択

bns_list = bns_scored_df.query('  score > 0.1 ').feature.unique().tolist()
corpus_data["bns_base_list"] = list(map(lambda text: " ".join(list(set(text.split(" ") ) & set(bns_list))) , corpus_data.text_wakati))

# word count feature
countvector = CountVectorizer()
matrix = countvector.fit_transform(corpus_data.bns_base_list)

X = matrix
y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

bns_list = bns_scored_df.query(' score > 0.1 ').feature.unique().tolist()

corpus_data["bns_base_list"] = list(map(lambda text: " ".join(list(set(text.split(" ") ) & set(bns_list))) , corpus_data.text_wakati))

# word count feature

countvector = CountVectorizer()

matrix = countvector.fit_transform(corpus_data.bns_base_list)

X = matrix

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

クラスタ、トピック、分散表現ベース

続いて、k-meansやLDAやword2vecを用いて特徴量を作成する方法です。今回はk-means、ミニバッチk-means、LDA、FastTextによる分散表現を扱います。

k-means、ミニバッチk-means

from sklearn.feature_extraction.text import CountVectorizer
from sklearn import preprocessing
from sklearn.cluster import  KMeans, MiniBatchKMeans

# word count feature
countvector = CountVectorizer(min_df=0.01, max_df=0.90)
matrix = countvector.fit_transform(corpus_data.text_wakati)

# 標準化
my_scalar = preprocessing.StandardScaler()
my_scalar.fit(pd.DataFrame(matrix.toarray()))
matrix_std = my_scalar.transform(pd.DataFrame(matrix.toarray()))

n_clusters = 15

kmeans = KMeans(n_clusters=n_clusters)
minibatch = MiniBatchKMeans(n_clusters=n_clusters)

kmeans.fit(matrix_std.transpose())
minibatch.fit(matrix_std.transpose())

# k-means
X = kmeans.cluster_centers_.transpose()
y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

# ミニバッチk-means
X = minibatch.cluster_centers_.transpose()
y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

from sklearn.feature_extraction.text import CountVectorizer

from sklearn import preprocessing

from sklearn.cluster import KMeans, MiniBatchKMeans

# word count feature

countvector = CountVectorizer(min_df=0.01, max_df=0.90)

matrix = countvector.fit_transform(corpus_data.text_wakati)

# 標準化

my_scalar = preprocessing.StandardScaler()

my_scalar.fit(pd.DataFrame(matrix.toarray()))

matrix_std = my_scalar.transform(pd.DataFrame(matrix.toarray()))

n_clusters = 15

kmeans = KMeans(n_clusters=n_clusters)

minibatch = MiniBatchKMeans(n_clusters=n_clusters)

kmeans.fit(matrix_std.transpose())

minibatch.fit(matrix_std.transpose())

# k-means

X = kmeans.cluster_centers_.transpose()

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

# ミニバッチk-means

X = minibatch.cluster_centers_.transpose()

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

LDA

こちらはgensimでLDAを推定し、推定したトピックの割合をデータフレームで返すコードです。

from gensim import corpora, models
import  gensim
import pickle

# データフレームのテキストデータをもとに、gensimで読み込み可能なフォーマットに変換するための関数
def nested_list(strings):
    words = strings.split(" ")
    return words

# トピック割合を格納するためのデータフレームを作成する関数
def making_topic_detaframe(integer, n):
    topic_table = pd.DataFrame(index=n.index)
    for topic_number in range(integer):
        column_name = "topic" + str(topic_number)
        topic_table[column_name] = 0
    return topic_table

# トピック割合を格納する関数
def topic_ratio_extract(corpus, topic_table):
    n_topic = len(topic_table.columns)
    i = 0
    for bow in corpus:
        t = lda.get_document_topics(bow)
        
        for each_topic in t:
            for topic_id in range(n_topic):
                if each_topic[0] == topic_id:
                    topic_table.iloc[i, topic_id]  = each_topic[1]
                
        i = i + 1
    return topic_table
 

corpus_data_filtered = corpus_data.query(' text_wakati !="" ').reset_index(drop=True)
corpus_data_filtered["text_wakati_fixed"] = list(map(lambda text:nested_list(text) ,corpus_data_filtered.text_wakati))

dictionary = gensim.corpora.Dictionary(corpus_data_filtered["text_wakati_fixed"])
dictionary.save_as_text('text.dict')

corpus = [dictionary.doc2bow(doc) for doc in corpus_data_filtered["text_wakati_fixed"]]
gensim.corpora.MmCorpus.serialize('text.mm', corpus)

tfidf = gensim.models.TfidfModel(corpus)
corpus_tfidf = tfidf[corpus]

with open('corpus_tfidf.dump', mode='wb') as f:
    pickle.dump(corpus_tfidf, f)

from gensim import corpora, models

import gensim

import pickle

# データフレームのテキストデータをもとに、gensimで読み込み可能なフォーマットに変換するための関数

def nested_list(strings):

words = strings.split(" ")

return words

# トピック割合を格納するためのデータフレームを作成する関数

def making_topic_detaframe(integer, n):

topic_table = pd.DataFrame(index=n.index)

for topic_number in range(integer):

column_name = "topic" + str(topic_number)

topic_table[column_name] = 0

return topic_table

# トピック割合を格納する関数

def topic_ratio_extract(corpus, topic_table):

n_topic = len(topic_table.columns)

i = 0

for bow in corpus:

t = lda.get_document_topics(bow)

for each_topic in t:

for topic_id in range(n_topic):

if each_topic[0] == topic_id:

topic_table.iloc[i, topic_id] = each_topic[1]

i = i + 1

return topic_table

corpus_data_filtered = corpus_data.query(' text_wakati !="" ').reset_index(drop=True)

corpus_data_filtered["text_wakati_fixed"] = list(map(lambda text:nested_list(text) ,corpus_data_filtered.text_wakati))

dictionary = gensim.corpora.Dictionary(corpus_data_filtered["text_wakati_fixed"])

dictionary.save_as_text('text.dict')

corpus = [dictionary.doc2bow(doc) for doc in corpus_data_filtered["text_wakati_fixed"]]

gensim.corpora.MmCorpus.serialize('text.mm', corpus)

tfidf = gensim.models.TfidfModel(corpus)

corpus_tfidf = tfidf[corpus]

with open('corpus_tfidf.dump', mode='wb') as f:

pickle.dump(corpus_tfidf, f)

トピック数をとりあえず30個くらいに指定して推定したトピックの割合を特徴量として文書分類を行います。そのため、特徴量の数は30個になります。

the_number_of_topic = 30

lda = gensim.models.ldamodel.LdaModel(
                            corpus=corpus,
                            alpha='asymmetric', 
                            eta = 'auto',
                            num_topics=the_number_of_topic,
                            id2word=dictionary,
                            random_state=123
)

# トピックの数に応じたデータフレームを作成
topic_table = making_topic_detaframe(the_number_of_topic, corpus_data_filtered)
# 推定したトピックごとの文書ごとの確率を格納
topic_table = topic_ratio_extract(corpus, topic_table)
# 元データとトピックテーブルを結合
corpus_data_filtered_with_topic = pd.concat([corpus_data_filtered, topic_table], axis=1)
# 実数に変換
corpus_data_filtered_with_topic[list(topic_table.columns)] = corpus_data_filtered_with_topic[list(topic_table.columns)].astype("float")


X = corpus_data_filtered_with_topic.filter(regex="topic").values
y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

the_number_of_topic = 30

lda = gensim.models.ldamodel.LdaModel(

corpus=corpus,

alpha='asymmetric',

eta = 'auto',

num_topics=the_number_of_topic,

id2word=dictionary,

random_state=123

)

# トピックの数に応じたデータフレームを作成

topic_table = making_topic_detaframe(the_number_of_topic, corpus_data_filtered)

# 推定したトピックごとの文書ごとの確率を格納

topic_table = topic_ratio_extract(corpus, topic_table)

# 元データとトピックテーブルを結合

corpus_data_filtered_with_topic = pd.concat([corpus_data_filtered, topic_table], axis=1)

# 実数に変換

corpus_data_filtered_with_topic[list(topic_table.columns)] = corpus_data_filtered_with_topic[list(topic_table.columns)].astype("float")

X = corpus_data_filtered_with_topic.filter(regex="topic").values

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

FastTextによる分散表現

今回はデータ数が少ないことから、学習済みの分散表現を用います。日本語のコーパスに対して、FastTextで推定された分散表現となります。学習済み分散表現はこちらから拝借しました。

import gensim.models.keyedvectors as word2vec_for_txt

# 少し時間がかかる
model_fasttext = word2vec_for_txt.KeyedVectors.load_word2vec_format('model.vec', binary=False)
w2v_fasttext = {w: vec for w, vec in zip(model_fasttext.wv.index2word, model_fasttext.wv.syn0)}

import gensim.models.keyedvectors as word2vec_for_txt

# 少し時間がかかる

model_fasttext = word2vec_for_txt.KeyedVectors.load_word2vec_format('model.vec', binary=False)

w2v_fasttext = {w: vec for w, vec in zip(model_fasttext.wv.index2word, model_fasttext.wv.syn0)}

分散表現は単語に対して計算されるので、単語に対して分散表現を足し合わせたものを特徴量として扱います。ここでは分散表現の合計値、平均値、TF-IDFで重みを付けた平均値の3つのパターンを試します。

合計値ベース

# 合計値ベース
X = [token.split(" ") for token in corpus_data.text_wakati]
X = np.array([ np.sum([w2v_fasttext[w] for w in words if w in w2v_fasttext] \
                        or [np.zeros(next(iter(w2v_fasttext.values())).size)], axis=0) for words in X])

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

# 合計値ベース

X = [token.split(" ") for token in corpus_data.text_wakati]

X = np.array([ np.sum([w2v_fasttext[w] for w in words if w in w2v_fasttext] \

or [np.zeros(next(iter(w2v_fasttext.values())).size)], axis=0) for words in X])

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

平均値ベース

# 平均値
X = [token.split(" ") for token in corpus_data.text_wakati]
X = np.array([ np.mean([w2v_fasttext[w] for w in words if w in w2v_fasttext] \
                        or [np.zeros(next(iter(w2v_fasttext.values())).size)], axis=0) for words in X])

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

# 平均値

X = [token.split(" ") for token in corpus_data.text_wakati]

X = np.array([ np.mean([w2v_fasttext[w] for w in words if w in w2v_fasttext] \

or [np.zeros(next(iter(w2v_fasttext.values())).size)], axis=0) for words in X])

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

TF-IDFで単語を重みづけた平均値ベース

# TF-IDFベース
from sklearn.feature_extraction.text import TfidfVectorizer
from collections import defaultdict

# TF-IDF値の計算
tfidf = TfidfVectorizer()
tfidf.fit(corpus_data.text_wakati)
max_idf = max(tfidf.idf_)

# TF-IDFの値を重み付けの係数として格納
word2weight = None
word2weight = defaultdict(lambda: max_idf, [(w, tfidf.idf_[i]) for w, i in tfidf.vocabulary_.items()])

# 分散表現の抽出
X = [token.split(" ") for token in corpus_data.text_wakati]
X = np.array([np.mean([w2v_fasttext[w] * word2weight[w] for w in words if w in w2v_fasttext] \
                      or [np.zeros(next(iter(w2v_fasttext.values())).size)], axis=0) for words in X ])

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape
print(n_features)

random_state = np.random.RandomState(0)
classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve
classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

# TF-IDFベース

from sklearn.feature_extraction.text import TfidfVectorizer

from collections import defaultdict

# TF-IDF値の計算

tfidf = TfidfVectorizer()

tfidf.fit(corpus_data.text_wakati)

max_idf = max(tfidf.idf_)

# TF-IDFの値を重み付けの係数として格納

word2weight = None

word2weight = defaultdict(lambda: max_idf, [(w, tfidf.idf_[i]) for w, i in tfidf.vocabulary_.items()])

# 分散表現の抽出

X = [token.split(" ") for token in corpus_data.text_wakati]

X = np.array([np.mean([w2v_fasttext[w] * word2weight[w] for w in words if w in w2v_fasttext] \

or [np.zeros(next(iter(w2v_fasttext.values())).size)], axis=0) for words in X ])

y = corpus_data.designer_flag.values

n_samples, n_features = X.shape

print(n_features)

random_state = np.random.RandomState(0)

classifier = RandomForestClassifier(random_state=random_state)

# classification and draw roc curve

classification_and_roc_analysis(k=10, classifier=classifier, X=X, y=y)

文書間の類似度ベース

今回は、デザイナーズマンションの定義文に似ているかどうかという観点で類似度ベースの特徴量を作ってみたいと思います。

from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.metrics.pairwise import cosine_similarity

designer_definition = "デザイナーズ物件というと、連想されがちなのが「コンクリート打ちっぱなしの外壁」に代表されるスタイリッシュなデザインの建物。\
                                        でも、それだけではありません。全く間仕切りがないワンルーム、螺旋階段が部屋の中にある部屋など、\
                                        普通の物件では考えられない設備やデザインが魅力です。また、キッチンが部屋の真ん中に配置されていたり、\
                                        バスやトイレがガラス張りであったり、洗濯機や電子レンジを収める収納スペースで生活感を出さない構造であったりすることも。\
                                        インテリアや家具までレイアウトされている部屋が多いのも特徴です。こうした一風変わった造りは、\
                                        デザイナーや建築家のこだわりやコンセプトが反映されています。そのため、デザイナーズ物件は建築家やデザインの\
                                        「作品」と言うこともできます。そうした作品に住むことは憧れであり、ある種のステイタスと言えるでしょう。"
designer_definition_wakati = nouns_extract(designer_definition) 


text_list = list(set(corpus_data.text_wakati))
text_list.extend(designer_definition_wakati )

# TF-IDFの計算
vectorizer = TfidfVectorizer()
vectorizer.fit_transform(text_list) # これで更新

# 類似度の計算
corpus_data["cos_sim"] = list(map(lambda number: cosine_similarity(vectorizer.transform([corpus_data.text_wakati[number]]),
                                                                                                                          vectorizer.transform([designer_definition_wakati]))[0][0] ,
                                                                                                                          range(corpus_data.index.size)))

bins = [0.00, 0.01, 0.03,0.05, 0.10,1.0]
labels=["0~1%", "1~3%", "3~5%", "5~10%", "10%~"]
corpus_data["cos_sim_bins"] = pd.cut(corpus_data.cos_sim, bins=bins, labels=labels, include_lowest=True)

df = corpus_data.groupby("cos_sim_bins")[["designer_flag"]].mean().reset_index()
ax = df.designer_flag.plot(xticks=df.index, rot=90, figsize=(10,5))
ax.set_xticklabels(df.cos_sim_bins);
ax.set_ylabel("designer_flag");
ax.set_xlabel("cos_sim_bins");

from sklearn.feature_extraction.text import TfidfVectorizer

from sklearn.metrics.pairwise import cosine_similarity

designer_definition = "デザイナーズ物件というと、連想されがちなのが「コンクリート打ちっぱなしの外壁」に代表されるスタイリッシュなデザインの建物。\

でも、それだけではありません。全く間仕切りがないワンルーム、螺旋階段が部屋の中にある部屋など、\

普通の物件では考えられない設備やデザインが魅力です。また、キッチンが部屋の真ん中に配置されていたり、\

バスやトイレがガラス張りであったり、洗濯機や電子レンジを収める収納スペースで生活感を出さない構造であったりすることも。\

インテリアや家具までレイアウトされている部屋が多いのも特徴です。こうした一風変わった造りは、\

デザイナーや建築家のこだわりやコンセプトが反映されています。そのため、デザイナーズ物件は建築家やデザインの\

「作品」と言うこともできます。そうした作品に住むことは憧れであり、ある種のステイタスと言えるでしょう。"

designer_definition_wakati = nouns_extract(designer_definition)

text_list = list(set(corpus_data.text_wakati))

text_list.extend(designer_definition_wakati )

# TF-IDFの計算

vectorizer = TfidfVectorizer()

vectorizer.fit_transform(text_list) # これで更新

# 類似度の計算

corpus_data["cos_sim"] = list(map(lambda number: cosine_similarity(vectorizer.transform([corpus_data.text_wakati[number]]),

vectorizer.transform([designer_definition_wakati]))[0][0] ,

range(corpus_data.index.size)))

bins = [0.00, 0.01, 0.03,0.05, 0.10,1.0]

labels=["0~1%", "1~3%", "3~5%", "5~10%", "10%~"]

corpus_data["cos_sim_bins"] = pd.cut(corpus_data.cos_sim, bins=bins, labels=labels, include_lowest=True)

df = corpus_data.groupby("cos_sim_bins")[["designer_flag"]].mean().reset_index()

ax = df.designer_flag.plot(xticks=df.index, rot=90, figsize=(10,5))

ax.set_xticklabels(df.cos_sim_bins);

ax.set_ylabel("designer_flag");

ax.set_xlabel("cos_sim_bins");

今回は変数が一つだけなので、機械学習はせず、デザイナーズマンション割合との関係を図示するにとどめておきます。横軸がデザイナーズマンションの定義と施設情報の類似度で、縦軸がデザイナーズマンション割合です。

どうやら、途中でデザイナーズマンション割合がピークを迎えるようです。

おわりに

最先端の手法は調べれていないですが、テキストデータを特徴量に落とし込む手段を備忘録として残しておきました。今回あげた中では、SOAベースの特徴量選択のAUCが83%と一番高かったですが、ベースラインが82%と僅差でした。そして、分散表現形のものは80%に届いた程度です。余力があれば新しい特徴量の作り方が分かり次第アップデートしようと思います。

追記

“Automate Stacking In Python How to Boost Your Performance While Saving Time”という記事を見つけたので、紹介されているvecstackモジュールを使って今回のモデルに関して簡単にstackingしてみようと思います。
コードに関しては、こちらのGitHubに上げています。試してみた所、AUCは88%になりました。結構上がりましたね。しかもコードはめちゃ短いので楽です。

参考文献

[1]Julian Avila et al(2019), 『Python機械学習ライブラリ scikit-learn活用レシピ80+』, impress top gear
[2]Receiver Operating Characteristic (ROC) with cross validation
[3]@Kensuke-Mitsuzawa(2016), “テキストデータで特徴量選択をする”, Qiita
[4]JapaneseTokenizer 1.6
[5]DocumentFeatureSelection 1.5
[6]自然言語処理における自己相互情報量 (Pointwise Mutual Information, PMI)
[7]【Techの道も一歩から】第3回「第11回テキストアナリティクス・シンポジウム」
[8]文書分類タスクでよく利用されるfeature selection

学習済み分散表現を用いた文書分類に挑戦（一部再学習も）

はじめに

2018年9月のテキストアナリティクスシンポジウムに行った際に、学習済みの分散表現で事前学習したモデルを使って分類してうまくいく事例が紹介されていました。
全てのタスクにおいてうまくいくとは思えませんが、試すコストはあまりかからないので試してみます。

2017年のテキストアナリティクスシンポジウムにおいても、メルカリやGunosyでは分散表現を用いた手法が一番精度が高いと言われていましたし、今年の会ではNLP系の学会でも分散表現はデファクトスタンダードになっているという話も伺いました。
2013~14年はLDAを使った研究が多かった気がしますが、徐々にシフトしていっているんですね。

これまで（Word2Vecを用いて蒙古タンメン中本の口コミ評価を予測してみる）は4000件程度の蒙古タンメン中本の口コミの情報を元に分散表現を手に入れていましたが、学習済みの分散表現を用いたアプローチも有効かもしれないと思い、試してみようと思います。

分類タスク

某グルメ口コミサイトの蒙古タンメン中本の口コミのテキストから、3.5点以上の評価かどうかを予測するタスクを扱います。
本当は、ポケモン図鑑の説明文から水やら炎やらのタイプを予測するとかをしたいのですが、あいにく手元にデータがないので、以前集めた蒙古タンメン中本の口コミを使います。（実は後日、ポケモン図鑑のデータを集めたのですが、平仮名にまみれたデータな上に、データ数も800件しかなかったので、どのみち厳しかったです。）

学習済み分散表現

Word2Vecなどで大量の文書をもとに学習させた分散表現のことを指します。
大規模コーパスで分散表現を手に入れる際は、数十GBにも相当するテキストデータを数時間かけて推定するので、学習済みのモデルは非常にありがたいです。（4年前に会社のPCで計算した際は、12時間くらいかかったこともありました。）

無料で提供してくださっている分散表現については、すでにこちらのブログで紹介されています。そこで紹介されているものに少し付け足すと、日本語の分散表現に関しては以下のようなものがあります。

白ヤギコーポレーションのモデル：Gensim
東北大学乾・岡崎研究室のモデル：Gensim
Facebookの学習済みFastTextモデル：Gensim
NWJC から取得した単語の分散表現データ (nwjc2vec)：Gensim
NNLM embedding trained on Google News：TensorFlow

そこで、今回は各種学習済み分散表現と蒙古タンメン中本コーパスで求めた分散表現の文書分類の性能バトルをしてみたいと思います。
ただ、分散表現ではなく、単語の頻度をもとに特徴量を作ったものが一番精度が高いのですが、分散表現同士の比較でもってどの学習済み分散表現が中本の口コミ分類に役に立ちそうなのかを明らかにしようと思います。（本来は分析という観点から即でボツですが、見苦しくも比較していきます。）

前処理

前処理は以下の通りで、テキストデータを分かち書きして、数値や低頻度・高頻度語を除外しています。

import pandas as pd

# データの読み込み
corpus_data = pd.read_pickle("nakamoto_corpus.pickle").reset_index(drop=True)

import collections
import MeCab
import mojimoji
from string import digits

remove_digits = str.maketrans('', '', digits)

tagger = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/")

def nouns_extract(line):
    keyword=[]
    node = tagger.parseToNode(line).next
    while node:
        if node.feature.split(",")[0] == "名詞":
            keyword.append(node.surface)
        node = node.next
    keyword = str(keyword).replace("', '"," ")
    keyword = keyword.replace("\'","")
    keyword = keyword.replace("[","")
    keyword = keyword.replace("]","")
    return keyword


#欠損データを除外する関数
def FilterNANData(dataset, column_name):
    result_drop_remove = copy.deepcopy(dataset)
    result_drop_remove = result_drop_remove[~result_drop_remove[column_name].isnull()].reset_index(drop=True)
    return result_drop_remove[column_name]

#形態素解析して名詞のみを抽出し単語の頻度を集計して降順で返す関数
def NounceSum(wordlist):
    m = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/")
    noun_list = [] # 重複を含めた名詞のリスト
    
    for i in wordlist:
        for l in m.parse (i).splitlines():
            try:
                if l != 'EOS' and l.split('\t')[1].split(',')[0] == '名詞': # EOSを除き名詞のみ抽出
                    noun_list.append(l.split('\t')[0]) # 見出し追加
            except:
                pass

    noun_cnt = collections.Counter(noun_list) # 各名詞の数え上げ
    nouns_data = pd.DataFrame.from_dict(noun_cnt, orient='index').reset_index()
    nouns_data.columns = ['nouns', 'count']
    nouns_data = nouns_data.sort_values(by=["count"], ascending=None)
    return nouns_data

#ストップワードを取り除く
def stop_word(documents):
    texts = [word for word in documents.lower().split() if word not in stoplist]
    texts = " ".join(texts)
    return texts

#リストでもらったテキストからストップワードを取り除く
def stop_word_tolist(documents):
    document_total = ' '.join(documents.tolist())
    texts = [[word for word in document_total.lower().split() if word not in stoplist]]
    return texts


#全角を半角にする
corpus_data["text_fixed"] =  list(map(lambda text: mojimoji.zen_to_han(text, kana=False) , corpus_data.text))
#数字を除外する
corpus_data["text_fixed"] =  list(map(lambda text: text.translate(remove_digits) , corpus_data.text_fixed))
#形態素解析する
corpus_data["text_wakati"] = list(map(lambda text:nouns_extract(text) , corpus_data.text_fixed))

#単語の頻度の計算
term_freq = pd.DataFrame(NounceSum(FilterNANData(corpus_data, "text_fixed")).reset_index(drop=True))
term_freq["ratio"] = term_freq["count"]/term_freq["count"].sum()

#不要語の除去（指定したものや、頻出のもの、頻度の低すぎるものを除外）
stoplist = ["ーー", "HP", "http://", "https://"]
stoplist2 = term_freq.query(' 10 > count | count > 1000  ').nouns.tolist()
stoplist.extend(stoplist2)
stoplist = set(stoplist)

corpus_data["text_wakati_fixed"] = list(map(lambda text:stop_word(text) , corpus_data.text_wakati))

#3.5点以上であれば1そうでなければ0
corpus_data["label"] = np.where(corpus_data.rating >= 3.5, 1, 0)

import pandas as pd

# データの読み込み

corpus_data = pd.read_pickle("nakamoto_corpus.pickle").reset_index(drop=True)

import collections

import MeCab

import mojimoji

from string import digits

remove_digits = str.maketrans('', '', digits)

tagger = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/")

def nouns_extract(line):

keyword=[]

node = tagger.parseToNode(line).next

while node:

if node.feature.split(",")[0] == "名詞":

keyword.append(node.surface)

node = node.next

keyword = str(keyword).replace("', '"," ")

keyword = keyword.replace("\'","")

keyword = keyword.replace("[","")

keyword = keyword.replace("]","")

return keyword

#欠損データを除外する関数

def FilterNANData(dataset, column_name):

result_drop_remove = copy.deepcopy(dataset)

result_drop_remove = result_drop_remove[~result_drop_remove[column_name].isnull()].reset_index(drop=True)

return result_drop_remove[column_name]

#形態素解析して名詞のみを抽出し単語の頻度を集計して降順で返す関数

def NounceSum(wordlist):

m = MeCab.Tagger("-d /usr/local/lib/mecab/dic/mecab-ipadic-neologd/")

noun_list = [] # 重複を含めた名詞のリスト

for i in wordlist:

for l in m.parse (i).splitlines():

try:

if l != 'EOS' and l.split('\t')[1].split(',')[0] == '名詞': # EOSを除き名詞のみ抽出

noun_list.append(l.split('\t')[0]) # 見出し追加

except:

pass

noun_cnt = collections.Counter(noun_list) # 各名詞の数え上げ

nouns_data = pd.DataFrame.from_dict(noun_cnt, orient='index').reset_index()

nouns_data.columns = ['nouns', 'count']

nouns_data = nouns_data.sort_values(by=["count"], ascending=None)

return nouns_data

#ストップワードを取り除く

def stop_word(documents):

texts = [word for word in documents.lower().split() if word not in stoplist]

texts = " ".join(texts)

return texts

#リストでもらったテキストからストップワードを取り除く

def stop_word_tolist(documents):

document_total = ' '.join(documents.tolist())

texts = [[word for word in document_total.lower().split() if word not in stoplist]]

return texts

#全角を半角にする

corpus_data["text_fixed"] = list(map(lambda text: mojimoji.zen_to_han(text, kana=False) , corpus_data.text))

#数字を除外する

corpus_data["text_fixed"] = list(map(lambda text: text.translate(remove_digits) , corpus_data.text_fixed))

#形態素解析する

corpus_data["text_wakati"] = list(map(lambda text:nouns_extract(text) , corpus_data.text_fixed))

#単語の頻度の計算

term_freq = pd.DataFrame(NounceSum(FilterNANData(corpus_data, "text_fixed")).reset_index(drop=True))

term_freq["ratio"] = term_freq["count"]/term_freq["count"].sum()

#不要語の除去（指定したものや、頻出のもの、頻度の低すぎるものを除外）

stoplist = ["ーー", "HP", "http://", "https://"]

stoplist2 = term_freq.query(' 10 > count | count > 1000 ').nouns.tolist()

stoplist.extend(stoplist2)

stoplist = set(stoplist)

corpus_data["text_wakati_fixed"] = list(map(lambda text:stop_word(text) , corpus_data.text_wakati))

#3.5点以上であれば1そうでなければ0

corpus_data["label"] = np.where(corpus_data.rating >= 3.5, 1, 0)

処理を施すとこのようなデータになります。

特徴量は、scikit-learnのCountVectorizerやTfidfVectorizer、分散表現の合計・平均・TF-IDFを求めたものを用意します。

from gensim.models.word2vec import Word2Vec
from gensim.models import word2vec
import matplotlib.pyplot as plt
import seaborn as sns

from tabulate import tabulate
from collections import Counter, defaultdict
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.pipeline import Pipeline
from sklearn.metrics import accuracy_score
from sklearn.cross_validation import cross_val_score

#単語の分散表現の合計値を求めるクラスの定義
class SumEmbeddingVectorizer(object):
    def __init__(self, word2vec):
        self.word2vec = word2vec
        self.dim = word2vec.values()
        self.dim = next(iter(self.dim))
        self.dim = self.dim.size
    
    def fit(self, X, y):
        return self 
 
    def transform(self, X):
        return np.array([
            np.sum([self.word2vec[w] for w in words if w in self.word2vec] 
                    or [np.zeros(self.dim)], axis=0)
            for words in X
        ])


#単語の分散表現の平均値を求めるクラスの定義
class MeanEmbeddingVectorizer(object):
    def __init__(self, word2vec):
        self.word2vec = word2vec
        self.dim = word2vec.values()
        self.dim = next(iter(self.dim))
        self.dim = self.dim.size
    
    def fit(self, X, y):
        return self 
 
    def transform(self, X):
        return np.array([
            np.mean([self.word2vec[w] for w in words if w in self.word2vec] 
                    or [np.zeros(self.dim)], axis=0)
            for words in X
        ])

#TF-IDFで重み付けした分散表現を求めるクラスの定義
class TfidfEmbeddingVectorizer(object):
    def __init__(self, word2vec):
        self.word2vec = word2vec
        self.word2weight = None
        self.dim = word2vec.values()
        self.dim = next(iter(self.dim))
        self.dim = self.dim.size
        
    def fit(self, X, y):
        tfidf = TfidfVectorizer(analyzer=lambda x: x)
        tfidf.fit(X)
        max_idf = max(tfidf.idf_)
        self.word2weight = defaultdict(
            lambda: max_idf, 
            [(w, tfidf.idf_[i]) for w, i in tfidf.vocabulary_.items()])
    
        return self
    
    def transform(self, X):
        return np.array([
                np.mean([self.word2vec[w] * self.word2weight[w]
                         for w in words if w in self.word2vec] or
                        [np.zeros(self.dim)], axis=0)
                for words in X
            ])

#入力変数と出力変数の指定
X, y = np.array(corpus_data.text_wakati_fixed), np.array(corpus_data.label)

#gensimで読み込むための形式にする
sentences = [token.split(" ") for token in corpus_data.text_wakati]

from gensim.models.word2vec import Word2Vec

from gensim.models import word2vec

import matplotlib.pyplot as plt

import seaborn as sns

from tabulate import tabulate

from collections import Counter, defaultdict

from sklearn.feature_extraction.text import CountVectorizer

from sklearn.feature_extraction.text import TfidfVectorizer

from sklearn.ensemble import ExtraTreesClassifier

from sklearn.pipeline import Pipeline

from sklearn.metrics import accuracy_score

from sklearn.cross_validation import cross_val_score

#単語の分散表現の合計値を求めるクラスの定義

class SumEmbeddingVectorizer(object):

def __init__(self, word2vec):

self.word2vec = word2vec

self.dim = word2vec.values()

self.dim = next(iter(self.dim))

self.dim = self.dim.size

def fit(self, X, y):

return self

def transform(self, X):

return np.array([

np.sum([self.word2vec[w] for w in words if w in self.word2vec]

or [np.zeros(self.dim)], axis=0)

for words in X

])

#単語の分散表現の平均値を求めるクラスの定義

class MeanEmbeddingVectorizer(object):

def __init__(self, word2vec):

self.word2vec = word2vec

self.dim = word2vec.values()

self.dim = next(iter(self.dim))

self.dim = self.dim.size

def fit(self, X, y):

return self

def transform(self, X):

return np.array([

np.mean([self.word2vec[w] for w in words if w in self.word2vec]

or [np.zeros(self.dim)], axis=0)

for words in X

])

#TF-IDFで重み付けした分散表現を求めるクラスの定義

class TfidfEmbeddingVectorizer(object):

def __init__(self, word2vec):

self.word2vec = word2vec

self.word2weight = None

self.dim = word2vec.values()

self.dim = next(iter(self.dim))

self.dim = self.dim.size

def fit(self, X, y):

tfidf = TfidfVectorizer(analyzer=lambda x: x)

tfidf.fit(X)

max_idf = max(tfidf.idf_)

self.word2weight = defaultdict(

lambda: max_idf,

[(w, tfidf.idf_[i]) for w, i in tfidf.vocabulary_.items()])

return self

def transform(self, X):

return np.array([

np.mean([self.word2vec[w] * self.word2weight[w]

for w in words if w in self.word2vec] or

[np.zeros(self.dim)], axis=0)

for words in X

])

#入力変数と出力変数の指定

X, y = np.array(corpus_data.text_wakati_fixed), np.array(corpus_data.label)

#gensimで読み込むための形式にする

sentences = [token.split(" ") for token in corpus_data.text_wakati]

蒙古タンメン中本の口コミ4000件から作成した分散表現：Gensim

まず、以前のブログで紹介した蒙古タンメン中本の分散表現ですが、以下のように推定しています。

#Word2Vecを実行する。 
model = Word2Vec(sentences,
                 sg=1,
                 size=50,
                 window=5, 
                 min_count=5,
                 workers=2,
                 seed=123)

#単語ごとの分散表現を手に入れる。
w2v = {w: vec for w, vec in zip(model.wv.index2word, model.wv.syn0)}

#Word2Vecを実行する。

model = Word2Vec(sentences,

sg=1,

size=50,

window=5,

min_count=5,

workers=2,

seed=123)

#単語ごとの分散表現を手に入れる。

w2v = {w: vec for w, vec in zip(model.wv.index2word, model.wv.syn0)}

Pipelineを用いてExtraTreesClassifierによる学習をします。特徴量は先程あげた、テキストベースのCountVectorizerやTfidfVectorizer、分散表現の合計・平均・TF-IDFで、評価指標はAUCのクロスバリデーションスコアとします。

#ベースラインとなる既存手法のモデルの準備
#etree
etree = Pipeline([("count_vectorizer",CountVectorizer(analyzer=lambda x: x)), 
                    ("extra trees", ExtraTreesClassifier(n_estimators=200))])
#etreeのTF-IDF版
etree_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)),
                          ("extra trees", ExtraTreesClassifier(n_estimators=200))])

#Word2Vecを特徴量としてExtraTreesによる分類器を準備する。
etree_w2v_sum = Pipeline([("word2vec vectorizer", SumEmbeddingVectorizer(w2v)), 
                        ("extra trees", ExtraTreesClassifier(n_estimators=200))])

etree_w2v = Pipeline([("word2vec vectorizer", MeanEmbeddingVectorizer(w2v)), 
                        ("extra trees", ExtraTreesClassifier(n_estimators=200))])
 
etree_w2v_tfidf = Pipeline([("word2vec vectorizer", TfidfEmbeddingVectorizer(w2v)), 
                        ("extra trees", ExtraTreesClassifier(n_estimators=200))])
 
#各モデルを実行し、クロスバリデーションスコアを計算し、出力させる。
all_models = [
    ("etree",etree),
    ("etree_tfidf",etree_tfidf),
    ("w2v_sum", etree_w2v_sum),    
    ("w2v", etree_w2v),
    ("w2v_tfidf", etree_w2v_tfidf)
]

scores = sorted([(name, cross_val_score(model, X, y, cv=5, scoring = "roc_auc").mean()) 
                 for name, model in all_models], 
                key = lambda x:x[0])

print(tabulate(scores, floatfmt=".4f", headers=("model", 'score')))

#ベースラインとなる既存手法のモデルの準備

#etree

etree = Pipeline([("count_vectorizer",CountVectorizer(analyzer=lambda x: x)),

("extra trees", ExtraTreesClassifier(n_estimators=200))])

#etreeのTF-IDF版

etree_tfidf = Pipeline([("tfidf_vectorizer", TfidfVectorizer(analyzer=lambda x: x)),

("extra trees", ExtraTreesClassifier(n_estimators=200))])

#Word2Vecを特徴量としてExtraTreesによる分類器を準備する。

etree_w2v_sum = Pipeline([("word2vec vectorizer", SumEmbeddingVectorizer(w2v)),

("extra trees", ExtraTreesClassifier(n_estimators=200))])

etree_w2v = Pipeline([("word2vec vectorizer", MeanEmbeddingVectorizer(w2v)),

("extra trees", ExtraTreesClassifier(n_estimators=200))])

etree_w2v_tfidf = Pipeline([("word2vec vectorizer", TfidfEmbeddingVectorizer(w2v)),

("extra trees", ExtraTreesClassifier(n_estimators=200))])

#各モデルを実行し、クロスバリデーションスコアを計算し、出力させる。

all_models = [

("etree",etree),

("etree_tfidf",etree_tfidf),

("w2v_sum", etree_w2v_sum),

("w2v", etree_w2v),

("w2v_tfidf", etree_w2v_tfidf)

]

scores = sorted([(name, cross_val_score(model, X, y, cv=5, scoring = "roc_auc").mean())

for name, model in all_models],

key = lambda x:x[0])

print(tabulate(scores, floatfmt=".4f", headers=("model", 'score')))

汗に関してコンテキストの似ている単語を抽出しています。

結果は、以下の通りで、分散表現を使わない方がAUCが高いです。ただ、w2v_tfidf（分散表現のTF-IDFを特徴量にしたもの）が分散表現の中でAUCが高いようです。今回はこの60.5%をベースラインに比較していこうと思います。

#中本の分散表現
model          score
-----------  -------
etree         0.6340
etree_tfidf   0.6499
w2v           0.5955
w2v_sum       0.5765
w2v_tfidf     0.6051

#中本の分散表現

model score

----------- -------

etree 0.6340

etree_tfidf 0.6499

w2v 0.5955

w2v_sum 0.5765

w2v_tfidf 0.6051

白ヤギコーポレーションのモデル：Gensim

こちらのリンク、「word2vecの学習済み日本語モデルを公開します」から、ダウンロードしてそのまま以下のコードでモデルを扱えます。

from gensim.models.word2vec import Word2Vec

model_path = "latest-ja-word2vec-gensim-model/word2vec.gensim.model"
model_shiroyagi = Word2Vec.load(model_path)
w2v_shiroyagi = {w: vec for w, vec in zip(model_shiroyagi.wv.index2word, model_shiroyagi.wv.syn0)}

from gensim.models.word2vec import Word2Vec

model_path = "latest-ja-word2vec-gensim-model/word2vec.gensim.model"

model_shiroyagi = Word2Vec.load(model_path)

w2v_shiroyagi = {w: vec for w, vec in zip(model_shiroyagi.wv.index2word, model_shiroyagi.wv.syn0)}

汗の関連語を抽出していますが、中国の歴史の何かですか？可汗とかいう単語は聞いたことあるかも。

まずは白ヤギさんの分散表現をそのまま使って予測してみます。（コードは先程のものとほぼ重複するので省略しています。）
残念ながら、ベースラインの60.5%には至りませんでした。

#学習済み分散表現をそのまま用いたもの
model          score
-----------  -------
w2v           0.6007
w2v_sum       0.5648
w2v_tfidf     0.5744

#学習済み分散表現をそのまま用いたもの

model score

----------- -------

w2v 0.6007

w2v_sum 0.5648

w2v_tfidf 0.5744

hogehoge.modelというフルモデル形式の場合は、再学習が可能です。詳しくはこちら（models.word2vec – Word2vec embeddings model）に書かれています。

train(sentences=None,
          corpus_file=None,
          total_examples=None, 
          total_words=None, 
          epochs=None, 
          start_alpha=None, 
          end_alpha=None, 
          word_count=0,
          queue_factor=2,
          report_delay=1.0,
          compute_loss=False,
          callbacks=())

train(sentences=None,

corpus_file=None,

total_examples=None,

total_words=None,

epochs=None,

start_alpha=None,

end_alpha=None,

word_count=0,

queue_factor=2,

report_delay=1.0,

compute_loss=False,

callbacks=())

今回は、白ヤギさんの分散表現に対して、追加で蒙古タンメン中本のテキストを食わせて再学習させます。

from gensim.models.word2vec import Word2Vec

model_path = "latest-ja-word2vec-gensim-model/word2vec.gensim.model"
model_shiroyagi = Word2Vec.load(model_path)

#リストの重複を除く
def remove_duplicates(l):
    return list(set(l))

#ユニーク単語のリストを作る
flat_list = [item for sublist in sentences for item in sublist]

model_shiroyagi.train(sentences=sentences,
                    total_examples=len(sentences),
                    total_words=len(remove_duplicates(flat_list)) ,
                    word_count= len(model_shiroyagi.wv.index2word),
                    epochs=4)

w2v_shiroyagi = {w: vec for w, vec in zip(model_shiroyagi.wv.index2word, model_shiroyagi.wv.syn0)}

from gensim.models.word2vec import Word2Vec

model_path = "latest-ja-word2vec-gensim-model/word2vec.gensim.model"

model_shiroyagi = Word2Vec.load(model_path)

#リストの重複を除く

def remove_duplicates(l):

return list(set(l))

#ユニーク単語のリストを作る

flat_list = [item for sublist in sentences for item in sublist]

model_shiroyagi.train(sentences=sentences,

total_examples=len(sentences),

total_words=len(remove_duplicates(flat_list)) ,

word_count= len(model_shiroyagi.wv.index2word),

epochs=4)

w2v_shiroyagi = {w: vec for w, vec in zip(model_shiroyagi.wv.index2word, model_shiroyagi.wv.syn0)}

ベースラインの60.5%よりも下回り、さきほどの白ヤギさんのもともとの分散表現よりも下回りました。

model          score
-----------  -------
w2v           0.5731
w2v_sum       0.5515
w2v_tfidf     0.5774

model score

----------- -------

w2v 0.5731

w2v_sum 0.5515

w2v_tfidf 0.5774

再学習してもかえって精度が下がったりすることから、簡単に精度が出るわけではなさそうです。まぁ、理想はその適用領域での大量のテキストデータがあることで、Wikipediaを元に作成した分散表現に強く依存しても駄目なのだろうと思われます。

東北大学乾・岡崎研究室のモデル：Gensim

日本語 Wikipedia エンティティベクトルからダウンロードした学習済み分散表現を用います。ダウンロード後は普通に「gzip -d file.txt.gz」みたいにターミナル上で解凍します。以下のコードを実行すればすぐに使うことができます。
ただし、KeyedVectors形式のものは白ヤギさんのように再学習ができません。（Why use KeyedVectors instead of a full model?）

import gensim.models.keyedvectors as word2vec_for_txt

#学習済み分散表現の読み込み
file_name = "jawiki.all_vectors.100d.txt"
model_tonpei = word2vec_for_txt.Word2VecKeyedVectors.load_word2vec_format(file_name)

＃単語の分散表現の取得
w2v_tonpei = {w: vec for w, vec in zip(model_tonpei.wv.index2word, model_tonpei.wv.syn0)}

import gensim.models.keyedvectors as word2vec_for_txt

#学習済み分散表現の読み込み

file_name = "jawiki.all_vectors.100d.txt"

model_tonpei = word2vec_for_txt.Word2VecKeyedVectors.load_word2vec_format(file_name)

＃単語の分散表現の取得

w2v_tonpei = {w: vec for w, vec in zip(model_tonpei.wv.index2word, model_tonpei.wv.syn0)}

汗の類似語に関しては、難しい単語が高めに出ているようです。

results = model_tonpei.most_similar(positive="汗", topn=10)

for result in results:
    print(result[0], '\t', result[1])

曷多   0.6981213092803955
拭い   0.6970207691192627
びっしょり    0.6847031116485596
垢    0.6723113059997559
喝    0.6681551933288574
裴羅   0.662636935710907
拭く   0.6584706902503967
洗う   0.6570826768875122
汚    0.6568220257759094
小便   0.6558053493499756

results = model_tonpei.most_similar(positive="汗", topn=10)

for result in results:

print(result[0], '\t', result[1])

曷多 0.6981213092803955

拭い 0.6970207691192627

びっしょり 0.6847031116485596

垢 0.6723113059997559

喝 0.6681551933288574

裴羅 0.662636935710907

拭く 0.6584706902503967

洗う 0.6570826768875122

汚 0.6568220257759094

小便 0.6558053493499756

残念ながら、ベースラインの60.5%には至りませんでした。

model

業界にいるものとしての思い

ジョブマッチング界隈の課題

ジョブマッチング界隈のデータの課題

ジョブマッチング系の論文色々

はじめに

目次

単純なRNNとは

LSTMとは

Bidirectional LSTMとは

Bidirectional LSTMで文書分類

データ

前処理

実行

追記

参考文献

はじめに

LDATSパッケージについて

データ

試してみた

おわりに

参考情報

はじめに

分析工程

データの収集

データの整形

単純集計・可視化

分析

APの発生件数の集計

人口あたりのAP件数

APの名前の集計

地価と人口あたりのAP発生件数の関係

APの詳細の集計

APの発生時期に関するテキストマイニング

APの詳細に関するテキストマイニング

おわりに

参考情報

テキストデータの特徴量化について

アプローチ

今回扱うデータ

今回扱うモデル

単語ベース

TF-IDFベースの特徴量選択

PMIベースの特徴量選択

SOAベースの特徴量選択

BNSベースの特徴量選択

クラスタ、トピック、分散表現ベース

k-means、ミニバッチk-means

LDA

FastTextによる分散表現

合計値ベース

平均値ベース

TF-IDFで単語を重みづけた平均値ベース

文書間の類似度ベース

おわりに

追記

参考文献

はじめに

分類タスク

学習済み分散表現

前処理

蒙古タンメン中本の口コミ4000件から作成した分散表現：Gensim

白ヤギコーポレーションのモデル：Gensim

東北大学 乾・岡崎研究室のモデル：Gensim

東北大学乾・岡崎研究室のモデル：Gensim