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が下がらないと過学習している可能性が高いので、これは過学習しているだけなのだろう。

参考文献

2019年に読んだデータ分析系の本の振り返り（21+1冊）

はじめに

2020年、あけましておめでとうございます。年末に自分自身を振り返ろうと思ったのですが、結局データ分析と勉強しかしていないわけで、書籍を振り返ろうと思うに至りました。私の知り合いのデータサイエンティストはだいたい全冊持っているであろうと思われますが、良い本だと思うので思い出していただければ幸いです。

1.『ベイズモデリングの世界』（岩波書店）

基本的に階層ベイズモデルを使って、個体ごとの異質性を考慮した分析手法が提案されています。前半はオムニバス形式で様々な先生がモデルの適用について執筆されており、後半では伊庭先生による階層ベイズモデルの講義になっています。途中でスタイン統計量による縮小推定の話があげられ、柔軟なモデリングのためには「階層化した方が少なくとも望ましい推定量が得られる」という数学的証明を捨てることもやむを得ないと書かれています。

2.『トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ) 』（コロナ社）

この本はトピックモデルの教科書というよりも、ベイズ推定の教科書という側面が強い印象があります。途中で出てくる数式は流し読みするのは難しく、最低2冊以上のノートが別途必要になると思います。一度でもLDAのパラメータを導出してみたいという方には良い教科書だと思います。疑似コードが提供されているので、それをもとにRやPythonでコーディングしていけば、一番シンプルなLDAが非常に短い行で実行できてしまうことに驚かれるかもしれません。人間が手を動かして推定アルゴリズムを導出しているからこそ、短いコードで済むということを実感できるはずです。

3.『構造的因果モデルの基礎』（共立出版）

グラフィカルなアプローチで因果推論を扱っている書籍です。Judea Pearl流の因果推論アプローチについて記すことを目的に書かれています。基礎と書かれていますが決して簡単ではありません。ただ、扱われる数学のレベルとしては確率と線形代数がわかれば大丈夫だと思われます。余談ではありますが、1章の相関関係と因果関係の事例紹介で「おむつとビールの話」が都市伝説ではなくきちんと記事としてWall Street Journalという雑誌に掲載されていたことが明らかにされています。

4.『現場で使える！PyTorch開発入門深層学習モデルの作成とアプリケーションへの実装 (AI & TECHNOLOGY)』（翔泳社）

PyTorchを触ったことがないが、深層学習の手法について知っている層を対象とした本です。6章まではGoogleのColabで動かせるのでGoogleに課金することなく深層学習による回帰、CNN、GAN、RNN、Encoder-Decoderモデル、ニューラル行列因子分解をPyTorchで試すことができます。写経したものはこちら。転移学習や高解像度化や画像生成、文章のクラス分類、文書生成、機械翻訳などもできるので、PyTorchでこれくらいの量をコーディングしたらこれくらいのことができるのかという学びや、他の人の書いたPyTorchコードを読みやすくなるなどの便益は十分にあると思いました。

5.『作ってわかる! アンサンブル学習アルゴリズム入門』（シーアンドアール研究所）

会社で行っているPythonもくもく会用に買った本で、scikit-learnを使わずに機械学習のアルゴリズム（アンサンブル系）をコーディングするための本です。pythonのコードについて逐次、細かい解説が行われているわけではないので、1行1行自分でコメントを加えながら写経をしていけば力が付くという本かなと思われます。sklearnはそれはそれで素晴らしいですが、こういう本でフルスクラッチで修行できるのはいいですね。

6.『数理統計学―基礎から学ぶデータ解析』（内田老鶴圃）

統計検定1級を合格された方のブログで紹介されていた教科書です。理系の大学生レベルの数学知識があれば、数理統計学の基礎を学べると思います。中心極限定理の証明や、様々な分布の期待値や分散、様々な分布の性質について数式を用いてしっかり理解することができます。数式もほどよく端折られているので、無論ですがノートが数冊必要になります。各章毎にある練習問題も解くことで力が付くと思います。日本の大学の授業の教科書がこれだったらジェノサイド（再履修者の大量発生）が起きるんだろうなと思ってしまった。

7.『44の例題で学ぶ統計的検定と推定の解き方』（オーム社）

統計の検定に関してだけ扱った珍しい本です。第3部までは統計学の普通の教科書ですが、それ以降であらゆる検定の例題が44件も載せられています。パラメトリックな検定から、ノンパラメトリックな検定まで幅広く扱われています。一番気にいっているのは仮説検定法の分類の表です。これさえあれば、どのデータに対してどの検定を行えばいいかが一目瞭然です。

8.『わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する』（技術評論社）

機械学習の原理を手早く数式を交えて学べる本です。かゆいところに手が届いていると言うか、既出の教科書では捨象されがちな、条件付き確率における2変数以上の条件づけでの表現に紙面を割いていたりしてくれるのが嬉しいです。ある程度数学の話はわかるが、だいぶ忘れているビジネスパーソンには大変にありがたいコンテンツと言えると思います。ベイズ線形回帰に関しても行列を用いた、わかりやすい導出方法が紹介されています。またコラムで紹介されている、測度論にどう向き合えばいいかの著者の見解は参考になります。

9.『Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)』

R言語とstanを用いてベイズ統計学を入門レベルから学べる本です。各トピックごとにそれなりの紙面が割かれています。例題も豊富にあるので、線形回帰・MCMC・情報量基準・階層ベイズモデルまで、ベイズ統計学を基礎から応用までしっかりと学べると思います。youtubeで著者の講義も配信されているので、留学気分を味わえます。

10.『scikit-learnとTensorFlowによる実践機械学習』（オライリージャパン）

2019年に日本で開かれたML SummitでTFの開発者がおすすめしていた教科書です。前半部分で機械学習の入門から応用までをわかりやすい説明で学ぶことができます。数式は少ないですが、図とソースコード（Python）がちりばめられており、手を動かして理解を進めることができます。後半部分はTensorFlowを用いた深層学習の基礎を同様に手を動かして学ぶことができます。ただ、TFのバージョンも変わってきているので前半の説明をアテにして読むのも良いと思います。

11.『AIアルゴリズムマーケティング自動化のための機械学習/経済モデル、ベストプラクティス、アーキテクチャ (impress top gear)』

マーケティングへのデータサイエンスの適用に関する珍しい書籍です。ソースコードはついていないですが、業務で使う際のアイデアが手に入ることもあります。一般的な回帰、生存時間分析、オークション、アトリビューション分析、アップリフトモデリング以外にも、情報検索やレコメンデーションやトピックモデルなどマーケティングながら学際的なトピックも扱われています。レコメンドなどで使われる、ランク学習に関して詳しく書かれた書籍をあまり知らないので、この本はその点においてもありがたい本でもあります。

12.『入門統計的因果推論』（朝倉書店）

ほぼ全ての章でグラフィカルなアプローチで因果推論を扱っています。例題も豊富なので、一つ一つ丁寧にやれば理解が捗ります。おそらく、例題の多さを含め一番丁寧にd分離性、do演算子、バックドア基準、フロントドア基準に関する説明をしてくれている本なのかなと思いました。グラフでの因果推論に関して初めての人でも、確率さえ知っていれば読み進めることができるはずです。また、途中で操作変数法の紹介もされ、経済学出身者としては読みやすい。ただ、傾向スコアのくだりや、DIDなどのくだりはあまり出てきません。あと、やってないですが章末の練習問題に対するSolution Manualが提供されているようです。

13.『実践ベイズモデリング -解析技法と認知モデル-』（朝倉書店）

ベイズモデリングを様々な事例に適用する方法がオムニバス形式で記された本です。ワイブル分布や異質性を考慮した二項分布、無制限複数選択形式のアンケートデータに対する手法、トピックモデル、項目反応理論などが扱われています。マーケティングの実務で使える事例が多いように感じました。こちらはサポートサイトでRコードとstanコードが提供されています。あと、appendixにあるプレート表現の見方も参考になります。

14.『機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)』

機械学習などで用いるベイズ推論を扱った教科書です。入門とありますが、入門者は書かれた数式をそのまま見ていても頭に入らないのではないでしょうか。手を動かしてなんぼの本だと思います。ノート2冊は絶対に必要です。たぶん、数式の展開を丁寧に記すと倍以上の厚みの本になると思います。各々のモデルに関してグラフィカルモデルが記されているのや、サンプルコードとしてGitHubにJuliaで書かれたソースコードが提供されているのも良いです。

15.『その問題、数理モデルが解決します』（ベレ出版）

物語形式で、様々な問題に対して数理モデリングのアプローチが紹介されています。途中でマッチング理論やゲーム理論やオークションなども登場することから、経済学出身者も喜ぶ内容かもしれません。社会人になってからナッシュ均衡という言葉が書かれた本は中々出会って来なかった。

16.『ヤバい予測学 ― 「何を買うか」から「いつ死ぬか」まであなたの行動はすべて読まれている』（CCCメディアハウス）

2013年と結構古い本ですが、データ分析を様々な事象に対して適用した事例紹介本です。アップリフトモデリングへの言及もあり、こういったものに関して日本は何年も遅れてブームが来るんだなという実感を与えてくれた本でもありました。appendixに分析事例が147個ほどあげられているのも参考になります。

17.『たのしいベイズモデリング2: 事例で拓く研究のフロンティア』（北大路書房）

主にstanを用いたベイズモデリングによる分析事例が1と2で38本もオムニバス形式で載っています。ほとんどの事例で階層ベイズモデルが扱われています。2では若干マーケティングに近い内容の題材も扱われ、データサイエンティストの人にも嬉しい内容かもしれません。もちろんデータとstanとRのコードがサポートサイトで提供されています。

18.『カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)』（共立出版）

状態空間モデルで時系列予測を行うための手法が記されている本です。RのKFASパッケージが全面に渡って扱われています。トレンドを考慮したり、カレンダー効果を追加したり、共変量を追加したりなど様々なアプローチが紹介されコードも伴っているわけですから、業務でも抜群に役に立ちました。

19.『機械学習のエッセンス -実装しながら学ぶPython,数学,アルゴリズム- (Machine Learning)』（SBクリエイティブ）

自分のいる会社で最低限の数学がわかると思われる若いメンバーに買ってもらうように言っている本です。微積分・線形代数だけでなく、カルシュ・キューン・タッカー条件（最適化数学）に関しても扱ってくれているので、ここで出てくる数学がわかれば大体の論文に立ち向かえると思います。さらに、Pythonの基礎もこれで学ぶことができるので一石二鳥な素敵な本ですね。また、最後の方でスクラッチでアルゴリズムを書くパートがあり、こちらも勉強になります。

20.『機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践 (オライリー・ジャパン)』(オライリー・ジャパン)

機械学習における前処理の指針を与えてくれる本です。Pythonのコードが提供されています。例えばですが、「テストデータにだけある、新しい単語は取り除いてしまえばいい」などの細かいアドバイスが何気に嬉しいです。「Effectコーディング」「特徴量ハッシング」「ビンカウンティング」「バックオフ」「leakage-proof統計量」などは読むまで知らないところだったので勉強になりました。

21.『データサイエンスのための統計学入門 ―予測、分類、統計モデリング、統計的機械学習とRプログラミング』（オライリージャパン）

データ分析の仕事をする上で最低限必要な知識を幅広く抑えることができる本です。数式は少なく、ところどころ出てくるコードはR言語です。参考文献などがブログだったりするため厳密さがめちゃあるわけではないですが、業務で使う分には問題ないと思います。分類問題において、AUCなどの評価指標だけでなく、予測値自体の探索的分析のすすめなどが書かれており、参考になりました。また、特徴量エンジンとしてのk-NN法の話も面白いと思いました。

[+α]『プログラマのためのGoogle Cloud Platform入門サービスの全体像からクラウドネイティブアプリケーション構築まで』（翔泳社）

Google Cloud Platformを初めて触るデータ分析者にはちょうど良い本です。説明もわかりやすいので、いきなりアカウントを作ってドキュメントを解読するよりかは戸惑いは減るはずです。この本を土台に、GCS・GCEを駆使してML系のAPIを呼び出して使うなどの最低限の操作は私でもできるようになりました。GCPの画面や機能もどんどん変わっていくので書籍を買ってもアレなんですが、歴史的な背景も若干記述されているので、それはそれで勉強になります。ただ、エンジニアにこの本を買うべきか聞いた際にネガティブな意見があったのですが、たぶん現役プログラマからすると簡単過ぎるからなんだろうなと思います。

終わりに

2019年もぼちぼち勉強できましたが、2020年もこれまで同様にノートとペンを大事にする勉強を続けていき、コーディングも分析ももっともっと数をこなして会社や社会に求められるようなデータ分析官を目指していこうと思います。あぁ、英会話などの勉強をする時間を作るのが難しい。

[Stan]ロジスティック回帰の階層ベイズモデルとk-foldsクロスバリデーション

はじめに

stanは意思決定のための分析などでのパラメータ推定に使うことが多く、機械学習のために扱うことはありませんでした。ただ、過去にリク面などでお話したデータサイエンティストの方はstanで機械学習していて、クロスバリデーションもしているとの発言をされていました。
先日、記事を漁っていたらstanでクロスバリデーションを行うためのコードを書いている方を見つけたので、その方のコードをもとに大人のirisであるwineデータを用いて、質の高いワインかどうかを分類するために階層ベイズでのロジスティック回帰モデルを回してみたいと思います。

データについて

UCI Machine Learning Repositoryにある、赤ワインの評価と成分のデータです。データに関する説明はワインの味（美味しさのグレード）は予測できるか？（１）で丁寧になされていますので、ご確認ください。今回は6点以上であれば1を、そうでなければ0を取るものを教師データとしています。

分析方針

今回は階層ベイズモデルを扱うことから、グループごとにロジスティック回帰のパラメータが異なるという仮定を置きます。そのために、citric.acidというデータを3つのカテゴリデータに変換して、それをグループとして扱います。モデルでは、今回のデータセットの変数を全て回帰項として使います。もちろん、回帰用の式からはcitric.acidは除外します。
全体の80%を訓練データに、20%をテストデータとして、10foldsクロスバリデーションでのstanによる予測結果の平均AUCを評価指標とします。最後に、テストデータを用いた予測のAUCを確かめます。また、階層ベイズモデルではないモデルでの10foldsクロスバリデーションでのAUCとも比較します

分析概要

・データ整形
・訓練データとテストデータの分割
・クロスバリデーション用のデータの作成
・stanの実行
・クロスバリデーション結果の出力
・テストデータでの予測
・非階層モデルとの比較

全体のコード以下のリンクにあります。
kick_logistic_regression_allowing_k_hold_cross_validation_hierachical.R
logistic_regression_allowing_k_fold_cross_validation_hierachical.stan

データ整形

階層ベイズで扱うグループをcitric.acidから作っています。

library(tidyverse)
library(rstan)
library(GGally)
library(shinystan)
library(pbmcapply)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

wine_dataset <- read.csv("dataset/winequality-red.csv")

# Visualization -----------------------------------------------------------
ggpairs(wine_dataset)


# Making group -----------------------------------------------------------
wine_dataset <- wine_dataset %>% mutate(citric_acid_group = if_else(citric.acid < 0.2, 1,
                                                                    if_else(citric.acid < 0.4, 2, 3)))

wine_dataset <- wine_dataset %>% select(-citric.acid)

library(tidyverse)

library(rstan)

library(GGally)

library(shinystan)

library(pbmcapply)

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

wine_dataset <- read.csv("dataset/winequality-red.csv")

# Visualization -----------------------------------------------------------

ggpairs(wine_dataset)

# Making group -----------------------------------------------------------

wine_dataset <- wine_dataset %>% mutate(citric_acid_group = if_else(citric.acid < 0.2, 1,

if_else(citric.acid < 0.4, 2, 3)))

wine_dataset <- wine_dataset %>% select(-citric.acid)

訓練データとテストデータの分割

ワインの質に関するバイナリーデータをこちらで作成し、80%を訓練データに、20%をテストデータに分割しています。

# Split Train and Test ----------------------------------------------------
smp_size <- floor(0.8 * nrow(wine_dataset))

## set the seed to make your partition reproducible
set.seed(123)
train_ind <- sample(seq_len(nrow(wine_dataset)), size = smp_size)

train <- wine_dataset[train_ind, ]
test <- wine_dataset[-train_ind, ]

y <- if_else(wine_dataset$quality > 5, 1, 0)
x <- as.matrix(wine_dataset %>% select(-quality))

y_train <- y[train_ind]
y_test <- y[-train_ind]

x_train <- x[train_ind,]
x_test <- x[-train_ind,]

x_train_group <- x_train[,ncol(x_train)]
x_test_group <- x_test[,ncol(x_test)]

x_train <- x_train[,1:(ncol(x_train)-1)]
x_test <- x_test[,1:(ncol(x_test)-1)]

x_train <- scale(x_train)
x_test <- scale(x_test)

# Split Train and Test ----------------------------------------------------

smp_size <- floor(0.8 * nrow(wine_dataset))

## set the seed to make your partition reproducible

set.seed(123)

train_ind <- sample(seq_len(nrow(wine_dataset)), size = smp_size)

train <- wine_dataset[train_ind, ]

test <- wine_dataset[-train_ind, ]

y <- if_else(wine_dataset$quality > 5, 1, 0)

x <- as.matrix(wine_dataset %>% select(-quality))

y_train <- y[train_ind]

y_test <- y[-train_ind]

x_train <- x[train_ind,]

x_test <- x[-train_ind,]

x_train_group <- x_train[,ncol(x_train)]

x_test_group <- x_test[,ncol(x_test)]

x_train <- x_train[,1:(ncol(x_train)-1)]

x_test <- x_test[,1:(ncol(x_test)-1)]

x_train <- scale(x_train)

x_test <- scale(x_test)

クロスバリデーション用のデータの作成

こちらのコードでは任意の数でクロスバリデーション用のデータを作成し、stanで扱う訓練用データのlistに追加しています。
また、参考にしているブログより転用したstan_kfoldという関数を定義しています。k分割した際のstanの推定結果をリストに格納するための関数です。

N <- length(y_train) # sample size
n_fold <- 10 # number of folds
K <- ncol(x_train) #n umber of predictors

# Create cross validation data --------------------------------------------
# create 10 folds of data
hh <- kfold_split_random(n_fold, N) #hh index the fold ID of each data point
holdout_k <- matrix(0, nrow = N, ncol = n_fold)
for(i in 1:N) holdout_k[i, hh[i]] <- 1

# turn into a list
holdout_k <- split(holdout_k,rep(1:ncol(holdout_k),each=nrow(holdout_k)))


# the basic data object
data_m <- list(N=N,
               K=K,
               M = length(unique(x_train_group)),
               citric_acid_group = x_train_group,
               X=x_train,
               y=y_train
               )

# create a list of data list
data_l <- rep(list(data_m),n_fold)
# add the holdout index to it
for(i in 1:n_fold) data_l[[i]]$holdout <- holdout_k[[i]]


# Define function ---------------------------------------------------------
# function to parrallelize all computations
# need at least two chains !!!
stan_kfold <- function(file, list_of_datas, chains, cores,...){
  library(pbmcapply)
  badRhat <- 1.1 # don't know why we need this?
  n_fold <- length(list_of_datas)
  model <- stan_model(file=file)
  # First parallelize all chains:
  sflist <- 
    pbmclapply(1:(n_fold*chains), mc.cores = cores, 
               function(i){
                 # Fold number:
                 k <- ceiling(i / chains)
                 s <- sampling(model, data = list_of_datas[[k]], 
                               chains = 1, chain_id = i)
                 return(s)
               })
  
  # Then merge the K * chains to create K stanfits:
  stanfit <- list()
  for(k in 1:n_fold){
    inchains <- (chains*k - (chains - 1)):(chains*k)
    #  Merge `chains` of each fold
    stanfit[[k]] <- sflist2stanfit(sflist[inchains])
  }  
  return(stanfit) 
}

N <- length(y_train) # sample size

n_fold <- 10 # number of folds

K <- ncol(x_train) #n umber of predictors

# Create cross validation data --------------------------------------------

# create 10 folds of data

hh <- kfold_split_random(n_fold, N) #hh index the fold ID of each data point

holdout_k <- matrix(0, nrow = N, ncol = n_fold)

for(i in 1:N) holdout_k[i, hh[i]] <- 1

# turn into a list

holdout_k <- split(holdout_k,rep(1:ncol(holdout_k),each=nrow(holdout_k)))

# the basic data object

data_m <- list(N=N,

K=K,

M = length(unique(x_train_group)),

citric_acid_group = x_train_group,

X=x_train,

y=y_train

)

# create a list of data list

data_l <- rep(list(data_m),n_fold)

# add the holdout index to it

for(i in 1:n_fold) data_l[[i]]$holdout <- holdout_k[[i]]

# Define function ---------------------------------------------------------

# function to parrallelize all computations

# need at least two chains !!!

stan_kfold <- function(file, list_of_datas, chains, cores,...){

library(pbmcapply)

badRhat <- 1.1 # don't know why we need this?

n_fold <- length(list_of_datas)

model <- stan_model(file=file)

# First parallelize all chains:

sflist <-

pbmclapply(1:(n_fold*chains), mc.cores = cores,

function(i){

# Fold number:

k <- ceiling(i / chains)

s <- sampling(model, data = list_of_datas[[k]],

chains = 1, chain_id = i)

return(s)

})

# Then merge the K * chains to create K stanfits:

stanfit <- list()

for(k in 1:n_fold){

inchains <- (chains*k - (chains - 1)):(chains*k)

# Merge `chains` of each fold

stanfit[[k]] <- sflist2stanfit(sflist[inchains])

}

return(stanfit)

}

stanの実行

こちらのstanのコードでは、M個のグループごとにパラメータが異なるというモデルを書いています。modelブロックの途中でholdoutを入れることで一部のデータを推定に使わないようにしています。

data {
  int<lower=0> N;   // number of data items
  int<lower=0> K;   // number of predictors
  int<lower=0> M;   // number of group
  row_vector[K] X[N];
  int<lower=0,upper=1> y[N];
  int<lower=0,upper=M> citric_acid_group[N];
  int<lower=0, upper=1> holdout[N];
  // index whether the observation should be held out (1) or used (0)
}

parameters {
  real mu[K];
  real<lower=0> sigma[K];
  vector[K] beta[M];
}

model {
  for (k in 1:K) {
    mu[k] ~ normal(0, 100);
    sigma[k] ~ inv_gamma(1, 1);

    for (m in 1:M)
      beta[m,k] ~ normal(mu[k], sigma[k]);
  }
  for (n in 1:N){
    if(holdout[n] == 0){
      target += bernoulli_lpmf( y[n] | inv_logit(X[n] * beta[citric_acid_group[n]]));
  }
 }
}

data {

int<lower=0> N; // number of data items

int<lower=0> K; // number of predictors

int<lower=0> M; // number of group

row_vector[K] X[N];

int<lower=0,upper=1> y[N];

int<lower=0,upper=M> citric_acid_group[N];

int<lower=0, upper=1> holdout[N];

// index whether the observation should be held out (1) or used (0)

}

parameters {

real mu[K];

real<lower=0> sigma[K];

vector[K] beta[M];

}

model {

for (k in 1:K) {

mu[k] ~ normal(0, 100);

sigma[k] ~ inv_gamma(1, 1);

for (m in 1:M)

beta[m,k] ~ normal(mu[k], sigma[k]);

}

for (n in 1:N){

if(holdout[n] == 0){

target += bernoulli_lpmf( y[n] | inv_logit(X[n] * beta[citric_acid_group[n]]));

}

こちらはstanをキックするためのコードです。いつもと違い、先程定義したstan_kfoldを用いています。

# Kick the stan code ------------------------------------------------------
# run the functions
ss <- stan_kfold(file="model/logistic_regression_allowing_k_fold_cross_validation_hierachical.stan",
                 data_l,
                 chains=4,
                 cores=2)

# Kick the stan code ------------------------------------------------------

# run the functions

ss <- stan_kfold(file="model/logistic_regression_allowing_k_fold_cross_validation_hierachical.stan",

data_l,

chains=4,

cores=2)

クロスバリデーション結果の出力

以下は、k個ずつ手に入ったクロスバリデーションでの推定結果から、各パラメータの平均値を計算し、ロジスティック回帰モデルで２値の予測を行い、平均AUCを計算するコードです。

# Calculate Mean AUC ------------------------------------------------------
## hierachical model
set.seed(123)
cv_mean_auc <- NULL
for (i in 1:n_fold){
  ext_fit <- extract(ss[[i]]) # choose 1 chunk
  
  # グループによって推定したパラメータが違う
  coef_list <- NULL
  group_list <- data_l[[i]]$citric_acid_group[data_l[[i]]$holdout > 0]
  for(j in 1:length(group_list)){
    coef_list <- rbind(coef_list, colMeans(ext_fit$beta[,group_list[j],]))
  }
  
  lin_comb <- rowSums(data_l[[i]]$X[data_l[[i]]$holdout > 0, ] * coef_list)
  prob <- 1/(1 + exp(-lin_comb))
  pred_value <- rbinom(sum(data_l[[i]]$holdout), 1, prob)
  
  # Syntax (response, predictor):
  auc = pROC::auc(data_l[[i]]$y[data_l[[i]]$holdout > 0], pred_value)[1]
  cv_mean_auc <- append(cv_mean_auc, auc)
}

cv_mean_auc
mean(cv_mean_auc)

# Calculate Mean AUC ------------------------------------------------------

## hierachical model

set.seed(123)

cv_mean_auc <- NULL

for (i in 1:n_fold){

ext_fit <- extract(ss[[i]]) # choose 1 chunk

# グループによって推定したパラメータが違う

coef_list <- NULL

group_list <- data_l[[i]]$citric_acid_group[data_l[[i]]$holdout > 0]

for(j in 1:length(group_list)){

coef_list <- rbind(coef_list, colMeans(ext_fit$beta[,group_list[j],]))

}

lin_comb <- rowSums(data_l[[i]]$X[data_l[[i]]$holdout > 0, ] * coef_list)

prob <- 1/(1 + exp(-lin_comb))

pred_value <- rbinom(sum(data_l[[i]]$holdout), 1, prob)

# Syntax (response, predictor):

auc = pROC::auc(data_l[[i]]$y[data_l[[i]]$holdout > 0], pred_value)[1]

cv_mean_auc <- append(cv_mean_auc, auc)

}

cv_mean_auc

mean(cv_mean_auc)

平均AUCは0.675となりました。すごくいいわけではないですが、手抜きモデルとしてはまずまずと言ったところでしょうか。

テストデータでの予測

以下のコードで最初に分けていたテストデータでの予測結果を返します。

# Prediction --------------------------------------------------------------
# Choose best model in cross-validation
ext_fit <- extract(ss[[1]]) # choose 1 chunk

# グループによって推定したパラメータが違う
beta_post <- NULL
group_list <- x_test_group

for(j in 1:length(group_list)){
  beta_post <- rbind(beta_post, colMeans(ext_fit$beta[,group_list[j],]))
}

lin_comb <- rowSums(x_test * beta_post)
prob <- 1/(1 + exp(-lin_comb))
pred_value <- rbinom(nrow(x_test), 1, prob)


# Syntax (response, predictor):
auc = pROC::auc(y_test, pred_value)[1]
auc

# Prediction --------------------------------------------------------------

# Choose best model in cross-validation

ext_fit <- extract(ss[[1]]) # choose 1 chunk

# グループによって推定したパラメータが違う

beta_post <- NULL

group_list <- x_test_group

for(j in 1:length(group_list)){

beta_post <- rbind(beta_post, colMeans(ext_fit$beta[,group_list[j],]))

}

lin_comb <- rowSums(x_test * beta_post)

prob <- 1/(1 + exp(-lin_comb))

pred_value <- rbinom(nrow(x_test), 1, prob)

# Syntax (response, predictor):

auc = pROC::auc(y_test, pred_value)[1]

auc

実行の結果、AUCは0.665と、クロスバリデーションでの平均AUCと比べてあまり下がりませんでした。

非階層モデルとの比較

非階層モデルでも同様に10foldsクロスバリデーションの平均AUCを計算しました。非階層モデルよりもAUCが1%ポイントくらいは高いようです。

> mean(cv_mean_auc)
[1] 0.6745282

> mean(cv_mean_auc_normal)
[1] 0.6640103

> mean(cv_mean_auc)

[1] 0.6745282

> mean(cv_mean_auc_normal)

[1] 0.6640103

おわりに

現時点において、stanでの柔軟なモデリングを機械学習に活かす作法について紹介されている文献はあまりなく、選手人口はどれくらいいるのか気になるところですが、発見したブログのやり方でクロスバリデーションをカジュアルに行えるので、より多くの方がstanでの機械学習にチャレンジしうるものだなと思いました。ただ、このレベルの階層ベイズだとrstanarmで簡単にできてしまうので、より深く分析してモデルをカスタムしていきたいですね。

参考情報

[1]Lionel Hertzog (2018), “K-fold cross-validation in Stan,datascienceplus.com”
[2]Alex Pavlakis (2018), “Making Predictions from Stan models in R”, Medium
[3]Richard McElreath (2016), “Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)”, Chapman and Hall/CRC
[4]松浦健太郎 (2016), 『StanとRでベイズ統計モデリング (Wonderful R)』, 共立出版
[5]馬場真哉 (2019), 『実践Data Scienceシリーズ RとStanではじめるベイズ統計モデリングによるデータ分析入門』, 講談社

[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

Uplift Modeling用のパッケージtools4upliftを使ってみた

はじめに

今回は、今後仕事で使いたいという思いもあり、RでUplift Modelingに関して便利なパッケージがないか探した結果、2019年に登場したばかりのtools4upliftの存在を知りました。アップリフトモデリングのモチベーションに関しても簡単に説明しながら、サンプルデータで実践してみようと思います。

・Uplift Modelingとはなにか
・Uplift Modelingの卑近な例え話
・Uplift Modelingのサンプルデータ
・tools4upliftについて
・tools4upliftでCriteoデータを試してみる
・『仕事ではじめる機械学習』の9章のコードをCriteoデータに試してみる
・おわりに
・参考文献

Uplift Modelingとはなにか

きちんとした説明は、あまりにも今更感があるので説明は端折りたいと思います。既出の文献がありますので、そちらを熟読ください。

Uplift Modelingの卑近な例え話

自分が吉野家のマーケティング担当だとしましょう。吉野家のアプリで割引クーポンを顧客にばらまくことができるとします。
マーケターとして重要なのは、割引クーポンを渡したことをきっかけとして吉野家に足を運び購入する顧客を増やせるかどうかになります。

マーケターの手元にあるのは、割引クーポンをばらまいた顧客とばらまかなかった顧客、そして吉野家で牛丼を食べたかどうかのデータです。
以前のマーケティング担当者がランダムにクーポンをばらまいていたことが重要なポイントです。

このデータから、顧客は以下の4分類に分かれます。

無関心：割引クーポンをばらまこうが我関せず。そもそも吉野家に行く気はない。
説得可能：普段、牛丼が安いすき屋にばかり行っているが、割高に感じている吉野家に負い目を感じている。割引クーポンで揺さぶられ来店する。
天の邪鬼：吉野家コピペのように、割引クーポンを握りしめた家族連れに遭遇したくないので、割引クーポンをばらまかれたら来店しないような客。
鉄板：毎日決まった時間に吉野家に行くことを心に決めている客。

マーケターは割引クーポンをばらまいた顧客と割引クーポンをばらまいていない顧客にデータを二分し、それぞれ機械学習のための訓練用データとテスト用データを用意します。

つまり、「割引クーポンをばらまいた顧客」の訓練用データとテスト用データと「割引クーポンをばらまいていない顧客」の訓練用データとテスト用データの計4つのデータセットを用意します。

まず、牛丼の購入の有無を教師とした訓練用データでロジスティック回帰モデルなどを推定します。
その結果、「割引クーポンをばらまいた顧客」から推定したモデルと、「割引クーポンをばらまいていない顧客」から推定したモデルが手元に残ります。

2つのテスト用データを1つにまとめて、先程推定したモデルを用いて、牛丼の購入確率を求めます。モデルは２つあるので、予測結果がテスト用データ1つに対して2つあることになります。

その予測結果の比（「割引クーポンをばらまいた顧客」モデルベースの予測値÷「割引クーポンをばらまいていない顧客」モデルベースの予測値）をアップリフトとみなします。

以下の図はこれまでの説明を図にしたものです。

アップリフトがどの程度の水準であれば、説得可能なユーザーが多いのかを探っていくことで、吉野家のアプリにおいて、どのユーザーに割引クーポンを発行するべきかがわかることになります。

Uplift Modelingのサンプルデータ

残念なことに吉野家のアプリのデータはありません。そこで今回は公開データを利用します。
以前より、The MineThatData E-Mail Analytics And Data Mining ChallengeのメールのデータがUplift Modelingで非常にしばしば取り上げられるデータでしたが、Twitterで他にデータないのかとぼやいたところ、2名の方にCriteo Uplift Prediction Datasetを紹介していただきました。

余談ですが、Criteo社と言えばディスプレイ広告のキング的な存在で、少し商品のリンクを踏んだだけであっという間に広告がレコメンドされますよね。自社で出稿用バナーを作っていましたが、CVRが高くなる良いクリエイティブを作ってきたのか、単にCriteo社のアルゴリズムが優秀なだけなのか非常に気になるところでしたね。

Criteo社が提供してくれている今回のデータは、2500万行に及ぶユーザーのデータで、プライバシー保護の観点から特徴量は復元できないような形式で提供されています。バイナリーのラベルとしては訪問やコンバージョンなどがあり、データ全体に占める処置群の割合は84.6%となっています。要は、吉野家で言う割引クーポンをばらまいた顧客が全体の84.6%に及ぶということです。

tools4upliftについて

2019年1月に公開されたRのUplift Modeling用のパッケージです。

特徴量における連続値をカテゴリ変数にする際に、最適な階級値を求めてくれる関数
アップリフトモデリングの可視化する関数
アップリフトモデリングにおける特徴量選択ができる関数
アップリフトモデリングにおけるモデルのバリデーションを行う関数

などが提供されており、ちょいとRを触れるマーケターにとって、アップリフトモデリングにおける試行錯誤がかなりしやすくなる便利なパッケージだと思いました。
なお、このパッケージで扱っているモデルはロジスティック回帰になります。介入データをもとに推定したモデルの条件付き確率と非介入データをもとに推定したモデルの条件付き確率の差をアップリフトとして推定しています。

このパッケージの解説論文においては、アップリフトモデリングの評価指標としてQini曲線というものが提案されていました。Qini曲線はローレンツ曲線のようなもので、Qini曲線とランダムに割り当てた際のアップリフト量の差分の合計をQini係数と定義しています。

tools4upliftでCriteoデータを試してみる

library(tools4uplift)
library(tidyverse)
library(data.table)

x = fread("gunzip -c criteo-uplift.csv.gz")

head(x,5)

# 1000万件でデータをサンプリング
sampleNum <- sample(nrow(x),10000000)
x_sub <- x[sampleNum,]

x_sub <- x_sub %>% select(-exposure,-visit)
summary(x_sub)

# Baseline models ---------------------------------------------------------
set.seed(123);
split.data1 <- SplitUplift(data = x_sub,
                           p = 0.7,
                           group = c("treatment", "conversion"))

train <- split.data1[[1]]
valid <- split.data1[[2]]

base.tm <- DualUplift(data = train,
                      treat = "treatment",
                      outcome = "conversion",
                      predictors = colnames(train[,1:12]))

# baseline model for control group
base.tm[[1]]

# baseline model for treatment group
base.tm[[2]]

# predict the uplift on the validation set
base.tm.valid <- DualPredict(data = valid,
                             treat = "treatment",
                             outcome = "conversion", 
                             model = base.tm, nb.group = 5)[[1]]

# evaluate the model performance
base.tm.perf <- QiniTable(data = base.tm.valid,
                          treat = "treatment",
                          outcome = "conversion",
                          prediction = "uplift_prediction",
                          nb.group = 5)

# Qini曲線の描画
QiniCurve(base.tm.perf, title = "")
# アップリフト量の棒グラフの描画
QiniBarPlot(base.tm.perf, title = "")
# Qini係数の算出
QiniArea(base.tm.perf)

library(tools4uplift)

library(tidyverse)

library(data.table)

x = fread("gunzip -c criteo-uplift.csv.gz")

head(x,5)

# 1000万件でデータをサンプリング

sampleNum <- sample(nrow(x),10000000)

x_sub <- x[sampleNum,]

x_sub <- x_sub %>% select(-exposure,-visit)

summary(x_sub)

# Baseline models ---------------------------------------------------------

set.seed(123);

split.data1 <- SplitUplift(data = x_sub,

p = 0.7,

group = c("treatment", "conversion"))

train <- split.data1[[1]]

valid <- split.data1[[2]]

base.tm <- DualUplift(data = train,

treat = "treatment",

outcome = "conversion",

predictors = colnames(train[,1:12]))

# baseline model for control group

base.tm[[1]]

# baseline model for treatment group

base.tm[[2]]

# predict the uplift on the validation set

base.tm.valid <- DualPredict(data = valid,

treat = "treatment",

outcome = "conversion",

model = base.tm, nb.group = 5)[[1]]

# evaluate the model performance

base.tm.perf <- QiniTable(data = base.tm.valid,

treat = "treatment",

outcome = "conversion",

prediction = "uplift_prediction",

nb.group = 5)

# Qini曲線の描画

QiniCurve(base.tm.perf, title = "")

# アップリフト量の棒グラフの描画

QiniBarPlot(base.tm.perf, title = "")

# Qini係数の算出

QiniArea(base.tm.perf)

こちらはアップリフト値の予測値の上位から右に並べた際のアップリフトの増大のグラフになります。20%あたりでピークになるようです。

こちらはアップリフト量の棒グラフです。20%の階級値を超えたらガクンと下がるのがわかります。

なお、Qini係数は0.03233551でした。

『仕事ではじめる機械学習』の9章のコードをCriteoデータに試してみる

tools4upliftの結果を鵜呑みにするのもあれなので、『仕事ではじめる機械学習』の9章のコードを使ってアップリフトモデリングを実践してみます。コードは丸パクリですが、謹んで掲載させていただきます。

%matplotlib inline

import pandas as pd
import matplotlib.pyplot as plt
from operator import  itemgetter
plt.style.use("ggplot")

from sklearn.model_selection import  train_test_split
from sklearn.linear_model import LogisticRegression

# https://ailab.criteo.com/criteo-uplift-prediction-dataset/
source_df = pd.read_csv("criteo-uplift.csv.gz")
source_df.head(10)

source_df.describe()

feature_vector_df = source_df.drop(["treatment","conversion","visit","exposure"],axis=1)

is_treat_list = list(source_df["treatment"] == 1)
is_cv_list = list(source_df["conversion"] == 1)

train_is_cv_list, test_is_cv_list, train_is_treat_list, \
        test_is_treat_list, train_feature_vector_df,\
        test_feature_vector_df = train_test_split(is_cv_list, is_treat_list,
                                                                                feature_vector_df,
                                                                                train_size=0.5,
                                                                                test_size=0.5,
                                                                                random_state=42)

treat_model = LogisticRegression(C=0.01)
control_model = LogisticRegression(C=0.01)

train_sample_num = len(train_is_cv_list)

treat_is_cv_list = [train_is_cv_list[i] for i in range(train_sample_num) if train_is_treat_list[i] == True]
treat_feature_vector_list = train_feature_vector_df[train_is_treat_list]

control_is_cv_list = [train_is_cv_list[i] for i in range(train_sample_num) if train_is_treat_list[i] == False]
control_feature_vector_list = train_feature_vector_df[list(map(lambda a:a == False ,train_is_treat_list))]

treat_model.fit(treat_feature_vector_list, treat_is_cv_list)
control_model.fit(control_feature_vector_list, control_is_cv_list)

treat_score = treat_model.predict_proba(test_feature_vector_df)
control_score = control_model.predict_proba(test_feature_vector_df)
score_list = treat_score[:,1] / control_score[:,1]

result = list(zip(test_is_cv_list, test_is_treat_list, score_list))
result.sort(key=itemgetter(2),reverse=True)
    

treat_uu = 0
control_uu = 0
treat_cv = 0
control_cv = 0
treat_cvr = 0.0
control_cvr = 0.0
lift = 0.0

stat_data = []

for is_cv, is_treat, score in result:
    if is_treat:
        treat_uu += 1
        if is_cv:
            treat_cv += 1
        treat_cvr = treat_cv / treat_uu
    else:
        control_uu += 1
        if is_cv:
            control_cv += 1
        control_cvr = control_cv / control_uu

    # コンバージョンレートの差に実験群の人数を掛けることでliftを算出
    lift = (treat_cvr - control_cvr) * treat_uu

    stat_data.append([is_cv, is_treat, score, treat_uu, control_uu, treat_cv, control_cv, treat_cvr, control_cvr, lift])


qdf = pd.DataFrame(columns=('treat_cvr', 'control_cvr'))

quantile_data = []
for n in range(10):
    start = int(n * len(result) / 10)
    end = int((n + 1) * len(result) / 10) - 1
    quantiled_result = result[start:end]
    
    treat_uu = list(map(lambda item:item[1], quantiled_result)).count(True)
    control_uu = list(map(lambda item:item[1], quantiled_result)).count(False)
    
    treat_cv = [item[0] for item in quantiled_result if item[1] == True].count(True)
    control_cv = [item[0] for item in quantiled_result if item[1] == False].count(True)
    
    treat_cvr = treat_cv / treat_uu
    control_cvr = control_cv / control_uu
    
    
    quantile_data.append([treat_uu, control_uu, treat_cv, control_cv, treat_cvr, control_cvr])
        
    label = "{}%~{}%".format(n*10, (n+1)*10)
    qdf.loc[label] = [treat_cvr, control_cvr]

qdf.plot.bar()
plt.xlabel("percentile")
plt.ylabel("conversion rate")

df = pd.DataFrame(stat_data)
df.columns = ["is_cv", "is_treat", "score", "treat_uu", "control_uu", "treat_cv", "control_cv", "treat_cvr", "control_cvr", "lift"]


# ベースラインを書き加える
df["base_line"] = df.index * df["lift"][len(df.index) - 1] / len(df.index)

df.plot(y=["treat_cv", "control_cv"])
plt.xlabel("uplift score rank")
plt.ylabel("conversion count")

df.plot(y=["treat_cvr", "control_cvr"], ylim=[0, 0.04])
plt.xlabel("uplift score rank")
plt.ylabel("conversion rate")

df.plot(y=["lift", "base_line"])
plt.xlabel("uplift score rank")
plt.ylabel("lift count")

df.plot(y=["treat_cv", "control_cv"], x="score", title="conversion count")
df.plot(y=["treat_cvr", "control_cvr"], ylim=[0, 0.04], x="score", title="conversion rate")
df.plot(y=["lift", "base_line"], x="score", title="lift")

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

%matplotlib inline

import pandas as pd

import matplotlib.pyplot as plt

from operator import itemgetter

plt.style.use("ggplot")

from sklearn.model_selection import train_test_split

from sklearn.linear_model import LogisticRegression

# https://ailab.criteo.com/criteo-uplift-prediction-dataset/

source_df = pd.read_csv("criteo-uplift.csv.gz")

source_df.head(10)

source_df.describe()

feature_vector_df = source_df.drop(["treatment","conversion","visit","exposure"],axis=1)

is_treat_list = list(source_df["treatment"] == 1)

is_cv_list = list(source_df["conversion"] == 1)

train_is_cv_list, test_is_cv_list, train_is_treat_list, \

test_is_treat_list, train_feature_vector_df,\

test_feature_vector_df = train_test_split(is_cv_list, is_treat_list,

feature_vector_df,

train_size=0.5,

test_size=0.5,

random_state=42)

treat_model = LogisticRegression(C=0.01)

control_model = LogisticRegression(C=0.01)

train_sample_num = len(train_is_cv_list)

treat_is_cv_list = [train_is_cv_list[i] for i in range(train_sample_num) if train_is_treat_list[i] == True]

treat_feature_vector_list = train_feature_vector_df[train_is_treat_list]

control_is_cv_list = [train_is_cv_list[i] for i in range(train_sample_num) if train_is_treat_list[i] == False]

control_feature_vector_list = train_feature_vector_df[list(map(lambda a:a == False ,train_is_treat_list))]

treat_model.fit(treat_feature_vector_list, treat_is_cv_list)

control_model.fit(control_feature_vector_list, control_is_cv_list)

treat_score = treat_model.predict_proba(test_feature_vector_df)

control_score = control_model.predict_proba(test_feature_vector_df)

score_list = treat_score[:,1] / control_score[:,1]

result = list(zip(test_is_cv_list, test_is_treat_list, score_list))

result.sort(key=itemgetter(2),reverse=True)

treat_uu = 0

control_uu = 0

treat_cv = 0

control_cv = 0

treat_cvr = 0.0

control_cvr = 0.0

lift = 0.0

stat_data = []

for is_cv, is_treat, score in result:

if is_treat:

treat_uu += 1

if is_cv:

treat_cv += 1

treat_cvr = treat_cv / treat_uu

else:

control_uu += 1

if is_cv:

control_cv += 1

control_cvr = control_cv / control_uu

# コンバージョンレートの差に実験群の人数を掛けることでliftを算出

lift = (treat_cvr - control_cvr) * treat_uu

stat_data.append([is_cv, is_treat, score, treat_uu, control_uu, treat_cv, control_cv, treat_cvr, control_cvr, lift])

qdf = pd.DataFrame(columns=('treat_cvr', 'control_cvr'))

quantile_data = []

for n in range(10):

start = int(n * len(result) / 10)

end = int((n + 1) * len(result) / 10) - 1

quantiled_result = result[start:end]

treat_uu = list(map(lambda item:item[1], quantiled_result)).count(True)

control_uu = list(map(lambda item:item[1], quantiled_result)).count(False)

treat_cv = [item[0] for item in quantiled_result if item[1] == True].count(True)

control_cv = [item[0] for item in quantiled_result if item[1] == False].count(True)

treat_cvr = treat_cv / treat_uu

control_cvr = control_cv / control_uu

quantile_data.append([treat_uu, control_uu, treat_cv, control_cv, treat_cvr, control_cvr])

label = "{}%~{}%".format(n*10, (n+1)*10)

qdf.loc[label] = [treat_cvr, control_cvr]

qdf.plot.bar()

plt.xlabel("percentile")

plt.ylabel("conversion rate")

df = pd.DataFrame(stat_data)

df.columns = ["is_cv", "is_treat", "score", "treat_uu", "control_uu", "treat_cv", "control_cv", "treat_cvr", "control_cvr", "lift"]

# ベースラインを書き加える

df["base_line"] = df.index * df["lift"][len(df.index) - 1] / len(df.index)

df.plot(y=["treat_cv", "control_cv"])

plt.xlabel("uplift score rank")

plt.ylabel("conversion count")

df.plot(y=["treat_cvr", "control_cvr"], ylim=[0, 0.04])

plt.xlabel("uplift score rank")

plt.ylabel("conversion rate")

df.plot(y=["lift", "base_line"])

plt.xlabel("uplift score rank")

plt.ylabel("lift count")

df.plot(y=["treat_cv", "control_cv"], x="score", title="conversion count")

df.plot(y=["treat_cvr", "control_cvr"], ylim=[0, 0.04], x="score", title="conversion rate")

df.plot(y=["lift", "base_line"], x="score", title="lift")

こちらの図はアップリフト値の階級値ごとのCVRです。最上位のアップリフト値はCVRの差が大きいですが、上位40~50%程度のアップリフト値のときにCVRの差が最も大きいようです。

アップリフト値の順位とCVRの図です。順位が低くても処置群のほうがCVRがわずかに高いようです。

アップリフトのスコアとCVRの関係です。2未満であればCVRは処置群が上回っていますが、一様な傾向はなさそうです。

コンバージョンレートの差に対象群の人数を掛けることでliftを算出したものです。アップリフトスコアが1~2点であれば儲かるようです。

tools4upliftと出している指標が違うので比較ができないのが難点に思いました。tools4upliftはオートマチックな感じで便利なのですが、『仕事ではじめる機械学習』の9章を正義として進めたいので、どうにか揃えれるようにしていきたいと思います。

おわりに

tools4upliftというマーケターにとって銀の弾丸になりそうなパッケージの存在を知ることができ、実際に非常に便利そうな関数が用意されているのがわかりました。ただ、開発されたばかりのパッケージなのでそこまで結果を信じていません。『仕事ではじめる機械学習』本の結果と揃えたいなと思いました。その点がはっきりすれば業務で使ってみるのも良いですし、任意のマーケターに安心して共有できると思います。

参考文献

[1] 有賀康顕・中山心太・西林孝 (2018) 『仕事ではじめる機械学習』オライリージャパン
[2] Mouloud Belbahri, Alejandro Murua, Olivier Gandouet, Vahid Partovi Nia (2019). “Uplift Regression: The R Package tools4uplift”, arXiv:1901.10867 [stat.AP]
[3] ohke (2019) 「Uplift modelingで施策が効く人を見極める」け日記
[4] usaito (2018) 「Uplift Modelingで介入効果を最適化する」 Qiita

ABEJA SIX 2019の1日目に行ってきましたレポート

今日は午後から有給をいただいて、ABEJA SIXの1日目に行ってきました。印象的だなと感じたものに関して、つらつらと雑記を載せておきたいと思います。

こちらは品川グランドプリンスホテルの庭園です。

こちらは会場の雰囲気です。

ブースの様子1です。

ブースの様子2です。

ABEJA SIX 2019

「食事画像認識モデル開発プロジェクトでの10個5個の教訓」株式会社FiNC Technologies 南野充則氏

FiNCは450万ダウンロードされているヘルスケア系のアプリを開発している会社。
ユーザーの継続率を高めるための施策として、機械学習を用いている。
今回の紹介事例ではユーザーの食事に関する情報を入力する手間を機械学習で短縮させ、短縮させることで継続率を高めることを狙っている。
食事の画像は1日に数万枚がアプリに投稿される。
食事の画像から栄養価などを計算することを目指している。
食事レシピ認識モデルでは、画像からレシピを識別し、メニューの量（グラム数）なども推定し、カテゴリ単位で決まっている栄養価から推定している。レシピ本の情報を入力したり、レシピサイトをクローリングし、レシピを一人あたりの栄養価になるように標準化などもしている。きれいな画像と栄養価（材料何グラムか）の伴ったクリーンなデータセットを用意するために自社のキッチンに料理人を呼び2000レシピ分の料理を作ったとのこと。
食材認識モデルでは食材一つ一つ（トマト一つとか、キャベツ一枚とか）を識別して、栄養価を素材単位で計算している。
学習の結果、管理栄養士よりも3%程度の誤差でメニューの栄養価を推定可能になった。
開発期間は6ヶ月間。
東大の松尾研にアドバイスをもらっているらしい。

5つの教訓

１．DL/ML人材をソフトウェアエンジニアから輩出すべき：
インフラ、サーバー、DB、パフォーマンスなどに明るいソフトウェアエンジニアが機械学習や深層学習を学ぶと、分析も実装もできる頼もしいメンバーになるので、ソフトウェアエンジニアのデータサイエンティスト化に注力しているらしい。目指すは論文のリプリケーションができるレベルとのこと。
2.データ取得から学習までのPDCAを最速にする：
ユーザーが画像を出したあとのフローをしっかりしていなかった。予期せぬデータが入ってくるので、そこへの対応も必要。アノテーションした項目を再学習するような仕組みを作り、そばの画像が苦手であれば、そばの画像を集中的に集めて学習させる。
3.オペレーションは自社で構築せよ：
泥臭い仕事と思い、丸投げしてはいけない。データセットの質が最も大事。データセットの質を担保するには評価手法を理解し細かいオペレーションを作る必要がある。アルバイトも自社で雇用、マネジャーもエンジニアとすることで当事者意識も芽生えやすい。
4.評価方法の決定からプロジェクトを始めよう：
AIを使えば、想像を超える何かが出てくると期待していまうフシがある。評価の仕方を決めたほうが、メンバーのゴールが見えるし。やりやすい。10%以内の誤差の難易度がどの程度なのかわからなかったりするし、解釈の多様性が生まれてしまうこともある。
5.プロジェクトはアジャイルで進めるべき：
作ったことのないモデルを作る際にスケジューリングを引くことは難しい。SOTAくらいいけますよと言ってしまい、自らを苦しめることになりかねない。

「機械学習におけるクラウド活用のポイント」アマゾンウェブサービスジャパン株式会社針原佳貴氏 & 宇都宮聖子氏

SageMakerいいぞというお話。
ビジネスにおいて、機械学習を進めるに際して重要なポイントは、
「ビジネス価値に落とし込む」
「データの流れを理解する」
「自分の力で頑張らない」
の3点が挙げられていた。
必要ではあるが、付加価値にはつながりにくい作業のことをUndifferentiated heavy liftingと呼ぶらしい。
機械学習プロジェクトを回す際に重要なこととして、
データ取得
↓
データ前処理
↓
モデルの開発・学習
↓
モデルの評価
↓
モデルの変換（エッジデバイスに送るにはデータを小さくする必要がある。）
↓
本番環境のデプロイ
↓
監視・評価データ変換
のループを繰り返すことが挙げられている。
S3(Simple Storage Service)に蓄積しているデータがあったとして、そのデータに対して、SageMakerで前処理やら機械学習を行い、学習済みの結果をS3にためれば、それを用いてエンドポイントの推論としてカジュアルに活用することができる。S3→SageMaker→S3のコンボが良いとのこと。
ここ1年間で200個くらいAWSのサービスやら機能が増えているので、それを知るだけでも大変そう。でもうまく使えば、Undifferentiated heavy liftingを避けることができる。
わからないことがあれば、ソリューションアーキテクトに質問したり、SageMakerのSlackで聞いたりすると良いらしい。
SageMakerでの学習の進め方としては3種類ある。１つ目は、TensorFlowなどでゴリゴリとアルゴリズムを書く。2つ目はAWS Marketplaceで販売されているアルゴリズムを時間単位で課金して使う。3つ目はAWSのビルトインのアルゴリズム（Object Detection、Semantic Segmentation、Factorization Machineなど）を使う。

「少数データからの学習法の展開とABEJAの取り組み」株式会社ABEJA 藤本敬介氏

データの質がモデルの結果を左右するが、きれいなデータを大量に集めるためにアノテーションをやるのは大変。少ないデータでも性能を出したい。
アプローチとしては、Data Augmentation、Transfer Learning、Meta learningの3つがある。

Data Augmentation（データ拡張）

データを擬似的に増やす手法
mixup:Beyond Empirical Risk Minimizationの研究では画像とラベルを合成してデータを水増しが行われている。
AutoAugment:LearningAugmentation Policies from Dataの研究では、強化学習で精度が上がりやすい水増し方法を学習している。
使わない手はない。

Transfer Learning（転移学習）

異なるデータセットで学習したものを再利用する。
Fine-tuning：別のデータで学習済みのモデルに対して、タスクに対してのデータに適用する。
Domain Adaptation：学習済みのモデルやデータの知識を再利用する。
Fine-tuningは有効な手段。

Meta learning

タスクの学習のしかたを学習する
少数のデータでのうまい学習方法を訓練しておいて、それを使い回す。

ABEJAの取り組み

データが少ない場合はFine-tuningで高精度を出しやすい。
External Network：中間層の情報を利用して、例外的な処理（ネットワークにバイパスみたいなものを通す）をすることで、Fine-tuningした際に精度が落ちないようにしている。不均衡データやクラス追加に対して強い手法とされている。データ数に応じてExternal Networkのサイズを調整でき、クラス1に大量のデータがある場合、1だけネットワークを深くして、2やら3はネットワークを浅くするなどの柔軟な対応が可能。これでもって不均衡データに対応できるとのこと。また、クラス追加に関しては、追加したクラスの分だけ学習すればいいようにネットワークの学習ができるらしい。ただし、学習に時間がかかるとのこと。
（よくわからないが）Model-Agnostic Meta-Learning（MAML（マムル））を応用したら精度が高まるらしい。

うーん、DNNは全然追いかけれていないので断片的にしかわからなかった。悔しいものです。

「Deep Learningの都市伝説と現実」株式会社ABEJA 白川達也氏

リサーチャーをする上で大事なこととしては、
1.先に見つけること
2.シンプルに解くこと
3.先に失敗する（大きな失敗は会社としてしないために）
の3つがある。
クリーンなデータで学習したほうが精度が高くなりやすく、過学習しにくい。ラベルの精度が高ければ、高いほどよい。Big Clean Data + DLで勝つる？
アノテーションは簡単ではない。アノテーターごとにわかりやすい情報がバラバラで、ブレるのが本質的。どこまでやるのか、どこが基準なのかというフレーミングとアンカーリングが重要。人間とかタスクを理解してすすめるのが良い。
半教師あり学習（アノテーションされていないデータを使って精度向上させる取り組み）も魅力的だが、教師データを増やしたほうが効率的。アノテーションできるならば、アノテーションしてしまおう。事前学習も意味があるので行う。
次にどんな技術がくるのか？ Graph Convolution、Annotation、Poincare Embeddings、ML in Hyperbolic Space
Taskonomyという研究が今後熱くなるかも。見たこともないタスクも解けるという柔軟性を持つモデルが構築できる？

感想

機械学習で精度を出すためにそこまで頑張るのか！という事例を聞けたり、知識として不足していたAWS系のサービスの話を聞けたり、自分の足りていない知識を補えた良いイベントだと思いました。

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

はじめに

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(mod

はじめに

目次

単純なRNNとは

LSTMとは

Bidirectional LSTMとは

Bidirectional LSTMで文書分類

データ

前処理

実行

追記

参考文献

はじめに

1.『ベイズモデリングの世界』（岩波書店）

2.『トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ) 』（コロナ社）

3.『構造的因果モデルの基礎』（共立出版）

4.『現場で使える！PyTorch開発入門 深層学習モデルの作成とアプリケーションへの実装 (AI & TECHNOLOGY)』（翔泳社）

5.『作ってわかる! アンサンブル学習アルゴリズム入門』（シーアンドアール研究所）

6.『数理統計学―基礎から学ぶデータ解析』（内田老鶴圃）

7.『44の例題で学ぶ統計的検定と推定の解き方』（オーム社）

8.『わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する』（技術評論社）

9.『Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)』

10.『scikit-learnとTensorFlowによる実践機械学習』（オライリージャパン）

11.『AIアルゴリズムマーケティング 自動化のための機械学習/経済モデル、ベストプラクティス、アーキテクチャ (impress top gear)』

12.『入門 統計的因果推論』（朝倉書店）

13.『実践 ベイズモデリング -解析技法と認知モデル-』（朝倉書店）

14.『機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)』

15.『その問題、数理モデルが解決します』（ベレ出版）

16.『ヤバい予測学 ― 「何を買うか」から「いつ死ぬか」まであなたの行動はすべて読まれている』（CCCメディアハウス）

17.『たのしいベイズモデリング2: 事例で拓く研究のフロンティア』（北大路書房）

18.『カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)』（共立出版）

19.『機械学習のエッセンス -実装しながら学ぶPython,数学,アルゴリズム- (Machine Learning)』（SBクリエイティブ）

20.『機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践 (オライリー・ジャパン)』(オライリー・ジャパン)

21.『データサイエンスのための統計学入門 ―予測、分類、統計モデリング、統計的機械学習とRプログラミング』（オライリージャパン）

[+α]『プログラマのためのGoogle Cloud Platform入門 サービスの全体像からクラウドネイティブアプリケーション構築まで』（翔泳社）

終わりに

はじめに

データについて

分析方針

分析概要

データ整形

訓練データとテストデータの分割

クロスバリデーション用のデータの作成

stanの実行

クロスバリデーション結果の出力

テストデータでの予測

非階層モデルとの比較

おわりに

参考情報

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

アプローチ

今回扱うデータ

今回扱うモデル

単語ベース

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

PMIベースの特徴量選択

SOAベースの特徴量選択

BNSベースの特徴量選択

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

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

LDA

FastTextによる分散表現

合計値ベース

平均値ベース

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

文書間の類似度ベース

おわりに

追記

参考文献

はじめに

Uplift Modelingとはなにか

Uplift Modelingの卑近な例え話

Uplift Modelingのサンプルデータ

tools4upliftについて

tools4upliftでCriteoデータを試してみる

『仕事ではじめる機械学習』の9章のコードをCriteoデータに試してみる

おわりに

参考文献

ABEJA SIX 2019

「食事画像認識モデル開発プロジェクトでの10個5個の教訓」 株式会社FiNC Technologies 南野 充則 氏

5つの教訓

4.『現場で使える！PyTorch開発入門深層学習モデルの作成とアプリケーションへの実装 (AI & TECHNOLOGY)』（翔泳社）

11.『AIアルゴリズムマーケティング自動化のための機械学習/経済モデル、ベストプラクティス、アーキテクチャ (impress top gear)』

12.『入門統計的因果推論』（朝倉書店）

13.『実践ベイズモデリング -解析技法と認知モデル-』（朝倉書店）

14.『機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)』

[+α]『プログラマのためのGoogle Cloud Platform入門サービスの全体像からクラウドネイティブアプリケーション構築まで』（翔泳社）

「食事画像認識モデル開発プロジェクトでの10個5個の教訓」株式会社FiNC Technologies 南野充則氏

「機械学習におけるクラウド活用のポイント」アマゾンウェブサービスジャパン株式会社針原佳貴氏 & 宇都宮聖子氏

「少数データからの学習法の展開とABEJAの取り組み」株式会社ABEJA 藤本敬介氏

「Deep Learningの都市伝説と現実」株式会社ABEJA 白川達也氏