概要

以前に投稿した記事
　Universal Sentence Encoderを使って文章の異常検知をする
では、Universal Sentence Encoder (USE)を用いて、夏目漱石の文章に混じった有価証券報告書の文章を見つけるというタスクを方向データの異常検知問題として扱いました。今回はUSEだけでなくELMoとBERTも用いて同種のタスクを解いて、３つのエンコーダーモデルを比較してみます。

日本語事前学習済みのELMoとBERTは、どちらもストックマークが公開しているモデルを使用します。

環境

計算はすべてGoogle Colaboratory上で行いました。BERTはTensorFlow 1.x系、USEはTensorFlow 2.x系を用いますので、環境を分けるため以下のような複数のノートブックで作業します。
data_preparation.ipynb ―― データ準備
ELMo_BERT_embedding.ipynb ―― ELMoとBERTによる埋め込みベクトルの計算
USE_embedding.ipynb ―― USEによる埋め込みベクトルの計算
anomaly_detection.ipynb ―― 異常検知の計算

データ準備

前回の記事では夏目漱石の小説の文章をメインにして有価証券報告書の文章を異常データとして混入させましたが、ストックマークのモデルはビジネスドメインのコーパスで事前学習されているので、今回は有価証券報告書をメインの文章とします。異常データはlivedoor newsコーパスから採取することとします。livedoor newsコーパスはnewsといっても通常のニュースではなく、家電やスポーツ・芸能ゴシップなどの記事から成るデータセットです。

まず、必要なライブラリをインポートして、グーグルドライブをマウントします。

data_preparation.ipynb

import re
import json
import glob
import numpy as np
from sklearn.model_selection import train_test_split

from google.colab import drive
drive.mount('/content/drive')

以下ではMy Drive下のanomaly_detectionというディレクトリで作業するものとします。適宜ディレクトリ名は読み替えてください。

data_preparation.ipynb

%cd /content/drive/'My Drive'/anomaly_detection

まず、chABSAデータセットをダウンロードし、有価証券報告書の文章のみを取り出します。手順は前回の記事と同様です。

data_preparation.ipynb

# データをダウンロードして展開
!wget https://s3-ap-northeast-1.amazonaws.com/dev.tech-sketch.jp/chakki/public/chABSA-dataset.zip
!unzip chABSA-dataset.zip
!rm chABSA-dataset.zip

# ファイルへのパスのリストを作成
chabsa_path_list = glob.glob("chABSA-dataset/*.json")

# 有価証券報告書の文章部分のみをリストに格納
chabsa_texts = []
for p in chabsa_path_list:
    with open(p, "br") as f:
        j =  json.load(f)
    for line in j["sentences"]:
        chabsa_texts += [line["sentence"].replace('\n', '')]

print(len(chabsa_texts))
# 6119

さらに短すぎる文章と長すぎる文章を削除します。

data_preparation.ipynb

def filter_by_length(texts_input, min_len=20, max_len=300):
    texts_output = []
    for t in texts_input:
        length = len(t)
        if length >= min_len and length <= max_len:
            texts_output.append(t)
    return texts_output

chabsa_texts = filter_by_length(chabsa_texts)
print(len(chabsa_texts))
# 5148

次にlivedoor newsコーパスをダウンロードします。何種類かある記事のうちsports-watchを使うこととします。

data_preparation.ipynb

# データをダウンロードして展開
!wget https://www.rondhuit.com/download/ldcc-20140209.tar.gz
!tar -xf ldcc-20140209.tar.gz
!rm ldcc-20140209.tar.gz

# ファイルへのパスのリストを作成
livedoor_path_list = glob.glob('./text/sports-watch/sports-watch-*.txt')
len(livedoor_path_list)
# 900

livedoor newsコーパスの文は■などの記号を多く含むので、リストアップした記号を除く前処理をします。また、URLを含む文は丸ごと除くことにします。さらに、コピーライトに関する行など英数字のみからなる行も除きます。

data_preparation.ipynb

def cleaning_text(texts):
    # URLを含む文を除く
    p = re.compile('https?://')
    if p.search(texts):
        return ''
    # 英数字のみの行も除く
    p = re.compile('[\w /:%#\$&\?\(\)~\.,=\+\-…（）＜＞]+')
    if p.fullmatch(texts):
        return ''
    # 改行や全角スペース、頻出の記号を除く
    remove_list = ['\n', '\u3000', '　', '＜', '＞', '・', 
                   '■', '□', '●', '○', '◆', '◇', 
                   '★', '☆', '▲', '△', '※', '＊', '*', '——']
    for s in remove_list:
        texts = texts.replace(s, '')
    return texts

livedoor_texts = []
for path in livedoor_path_list:
    with open(path) as f:
        texts = f.readlines()
    livedoor_texts += [cleaning_text(t) for t in texts[4:]] # 最初の３行は不要
print(len(livedoor_texts))

chABSAデータの文と形式を合わせるために句点を区切りに文章を分割し、文章の長さでフィルターをかけます。

data_preparation.ipynb

def split_texts(texts_input, split_by='。'):
    texts_output = []
    for t in texts_input:
        texts_output += t.split(split_by)
    return texts_output

livedoor_texts = split_texts(livedoor_texts)
print(len(livedoor_texts))
# 17522

livedoor_texts = filter_by_length(livedoor_texts)
print(len(livedoor_texts))
# 8149

２種類の文章群のリストが得られましたので、モデル開発用データとテスト用データを作成します。有価証券報告書の文章にラベル０を、livedoor newsの文章にラベル１を付与します。有価証券報告書の文章のうち８割を開発データに使用し、残り２割はテストデータに回します。開発データには約1%だけlivedoor newsの文章を混ぜます。テストデータは精度の評価がしやすいように２種類の文章の量は50%ずつとします。

data_preparation.ipynb

def make_dataset(main_texts, anom_texts, main_dev_rate=0.8, anom_dev_rate=0.01):
    len1 = len(main_texts)
    len2 = len(anom_texts)
    num_dev1 = int(len1 * main_dev_rate)
    num_dev2 = int(num_dev1 * anom_dev_rate)
    num_test1 = len1 - num_dev1
    num_test2 = num_test1

    print("開発データ　 main: {}, anom: {}".format(num_dev1, num_dev2))
    print("テストデータ main: {}, anom: {}".format(num_test1, num_test2))

    main_arr = np.hstack([np.reshape(np.zeros(len1), (-1, 1)), 
                          np.reshape(main_texts, (-1, 1))])
    anom_arr = np.hstack([np.reshape(np.ones(len2), (-1, 1)), 
                          np.reshape(anom_texts, (-1, 1))])
    
    dev1, test1 = train_test_split(main_arr, train_size=num_dev1)
    dev2, test2 = train_test_split(anom_arr, train_size=num_dev2)
    test2, _ = train_test_split(test2, train_size=num_test2)

    dev_arr = np.vstack([dev1, dev2])
    np.random.shuffle(dev_arr)
    test_arr = np.vstack([test1, test2])
    np.random.shuffle(test_arr)
    return dev_arr, test_arr

dev_arr, test_arr = make_dataset(chabsa_texts, livedoor_texts)
# 開発データ　 main: 4118, anom: 41
# テストデータ main: 1030, anom: 1030

print(dev_arr.shape, test_arr.shape)
# (4159, 2) (2060, 2)

作成したデータセットを保存します。

data_preparation.ipynb

np.save('dev.npy', dev_arr)
np.save('test.npy', test_arr)

テストデータをいくつかサンプルしてみると以下のようになります。

1: 仮にそのまま西武に逆転負けを喫していたなら、間違いなく大炎上となっていたであろう楽天の怠慢プレー
0: 資金運用収支は、有価証券利息配当金が減少したものの、貸出金利息が増加したことより、179億円（同8.4％増）となりました
1: 18歳の時から知ってますけど、子供だった彼が今や本当に大人になってきた
0: 事業承継ナビゲーターの設立により、比較的に長いスパンでいくつかの選択肢の中からベストな事業承継プランを選択したいというミッドキャップ（中堅企業）の経営者の方々のニーズに応えるべく営業活動を行っております
0: また、セグメント利益につきましては、エアバッグインフレータに起因する品質関連費用および米国の金利上昇に伴う販売費を中心とした諸経費等の増加、為替変動の影響、試験研究費の増加により、3,977億円と前連結会計年度に比べ1,460億円（26.8％）の減益となりました
1: それと、俺らやってきたんだから一回は監督に認めて貰いたいよなっていうのがあったんじゃないかな
1: そして、1位は「完璧にイメージ通りに蹴れたのは、あれだけかな
0: 航空輸入貨物の取扱いは堅調に推移したものの、為替の影響等により、売上高は1,017億円と前連結会計年度に比べ133億円、11.6％の減収となり、営業利益は11億円と前連結会計年度に比べ５億円、33.5％の減益となりました
0: 日本経済におきましては、雇用や所得環境の改善が続く中、景気には一部に改善の遅れも見られるものの緩やかな回復基調が継続しました
0: ＰＣコンテンツ配信事業につきましては、アーティスト及びタレント等のＰＣ向け有料ファンクラブサイトの運営、並びにオフィシャルサイトの受託制作などを実施し、他の事業部門も含め、将来の新たな収益獲得へと繋がることを見据えた事業展開を行ってまいりました

MeCab+NEologd辞書インストール

ここからノートブックELMo_BERT_embedding.ipynbで作業します。はじめに先ほど同様にGoogle Driveをマウントしておきます。

ELMo_BERT_embedding.ipynb

from google.colab import drive
drive.mount('/content/drive')

ストックマークの事前学習モデルはトークナイザーにMeCab+NEologd辞書が使用されています。MeCab+NEologd辞書のインストール方法は公式の配布ページに従います。まず、ノートブック上で以下のコマンドを実行してMeCab本体と他の必要なパッケージをインストールします。

ELMo_BERT_embedding.ipynb

!apt install swig mecab libmecab-dev mecab-ipadic-utf8 git make curl xz-utils file -y

次にNEologd辞書をインストールします。ディレクトリ'My Drive'下で作業するとディレクトリ名がスペースを含むためエラーが出るのでルートディレクトリで作業します。

ELMo_BERT_embedding.ipynb

%cd /content
!git clone --depth 1 https://github.com/neologd/mecab-ipadic-neologd.git
%cd mecab-ipadic-neologd
!echo yes | ./bin/install-mecab-ipadic-neologd -n
%cd /content/drive/'My Drive'/anomaly_detection2

最後にpythonでMeCabを呼び出すためのライブラリをインストールします。

ELMo_BERT_embedding.ipynb

!pip install mecab-python3
import MeCab

（2020/7/3 追記）記事の投稿時にはこれでうまく行っていたのですが、いつの間にかmecab-python3の実行時にエラーが出るようになってしまいました。エラーメッセージに従って、unidic-liteをインストールするとmecab-python3を実行できるようになります。

ELMo_BERT_embedding.ipynb

!pip install unidic-lite

今回使う日本語版ELMoではMeCabにより分かち書きしたトークンを入力する必要がありますので、先ほど作成したデータセットを読み込んでトークン化します。

ELMo_BERT_embedding.ipynb

%cd /content/drive/'My Drive'/anomaly_detection
dev_arr = np.load('dev.npy')
test_arr = np.load('test.npy')

# 文章のみ切り出す
dev_texts = dev_arr[:, 1]
test_texts = test_arr[:, 1]

def MeCab_tokenizer(texts):
    mecab = MeCab.Tagger(
        "-Owakati -d /usr/lib/x86_64-linux-gnu/mecab/dic/mecab-ipadic-neologd")
    token_list = []
    for s in texts:
        parsed = mecab.parse(s).replace('\n', '')
        token_list += [parsed.split()]
    return token_list

dev_tokens = MeCab_tokenizer(dev_texts)
test_tokens = MeCab_tokenizer(test_texts)

また、今回使用するBERTのプログラムではテキストファイルから入力する必要があるので、ファイルに書き込んでおきます。こちらは事前にトークン化する必要はありません。

ELMo_BERT_embedding.ipynb

with open('dev_text.txt', mode='w') as f:
    f.write('\n'.join(dev_texts))
with open('test_text.txt', mode='w') as f:
    f.write('\n'.join(test_texts))

ELMo

ELMoモデルはこちらの実装をストックマークの事前学習済みパラメータとともに使用します。具体的な使用方法は

大規模日本語ビジネスニュースコーパスを学習したELMo（MeCab利用）モデルの利用方法と精度比較検証

を参照しました。
まず、必要なライブラリをインストールし、リポジトリをクローンします。

ELMo_BERT_embedding.ipynb

!pip install overrides
!git clone https://github.com/HIT-SCIR/ELMoForManyLangs.git

setup.pyを実行することでインストールが完了します。

ELMo_BERT_embedding.ipynb

%cd ./ELMoForManyLangs
!python setup.py install
%cd ..

次に、こちらのリンクからストックマークの事前学習済みモデルをダウンロードします。モデルは「単語単位埋め込みモデル」と「文字単位・単語単位埋め込みモデル」があります。今回は「単語単位埋め込みモデル」を使用することにします。ダウンロードしたファイルをマウントしたGoogle Driveのフォルダ./ELMo_ja_word_levelに配置しました。フォルダの中身は以下のようになっているはずです。

ELMo_BERT_embedding.ipynb

!ls ./ELMo_ja_word_level/
# char.dic  config.json  configs  encoder.pkl  token_embedder.pkl  word.dic

これでモデルの準備が整いました。Embedderインスタンスを作成し、文章の埋め込みベクトルを計算するための関数を定義します。

ELMo_BERT_embedding.ipynb

from ELMoForManyLangs.elmoformanylangs import Embedder
from overrides import overrides

elmo_model_path = "./ELMo_ja_word_level"
elmo_embedder = Embedder(elmo_model_path, batch_size=64)

def get_elmo_embeddings(token_list, batch_size=64):
    length = len(token_list)
    n_loop = int(length / batch_size) + 1
    sent_emb = []
    for i in range(n_loop):
        token_emb = elmo_embedder.sents2elmo(
            token_list[batch_size*i: min(batch_size*(i+1), length)])
        for emb in token_emb:
          # sum over tokens to obtain the embedding for a sentence
          sent_emb.append(sum(emb))
    return np.array(sent_emb)

モデルの出力は入力トークンそれぞれに対する埋め込みベクトルです。文章内のトークンのベクトルすべての和を取ったものを文章の埋め込みベクトルとしています。ELMoベクトルの次元は1024です。

ELMo_BERT_embedding.ipynb

dev_elmo_embeddings = get_elmo_embeddings(dev_tokens)
test_elmo_embeddings = get_elmo_embeddings(test_tokens)
print(dev_elmo_embeddings.shape, test_elmo_embeddings.shape)
# (4159, 1024) (2060, 1024)

np.save('dev_elmo_embeddings.npy', dev_elmo_embeddings)
np.save('test_elmo_embeddings.npy', test_elmo_embeddings)

BERT

BERTに関してもストックマークが公開している事前学習済みモデルを使用します。こちらのダウンロードリンクからTensorFlow版をダウンロードし、フォルダ./BERT_base_stockmarkに配置しました。フォルダの中身は以下の通りです。

ELMo_BERT_embedding.ipynb

!ls ./BERT_base_stockmark
# bert_config.json  vocab.txt  output_model.ckpt.index  output_model.ckpt.meta output_model.ckpt.data-00000-of-00001

バージョン1.x系を指定してtensorflowをインポートします。

ELMo_BERT_embedding.ipynb

%tensorflow_version 1.x
import tensorflow as tf

公式リポジトリからモデル本体をクローンします。

ELMo_BERT_embedding.ipynb

!git clone https://github.com/google-research/bert.git

公式リポジトリには埋め込みベクトルを取り出すためのスクリプトコードextract_features.pyが用意されているので、オリジナルの英語版ならばそれをそのまま実行すればいいのですが、日本語事前学習モデルを使うためにはいくつかのファイルを編集する必要があります。

こちらのページの利用方法の説明に従って、tokenization.pyをトークナイザーにMeCabを使用するように変更します。まず、未知語に対して[UNK]のidである1を返すように、関数convert_tokens_to_ids(vocab, tokens)を以下のように変更します。

bert/tokenization.py

def convert_tokens_to_ids(vocab, tokens):
  # modify so that it returns id=1 which means [UNK] when a token is not in vocab
  output = []
  for t in tokens:
    if t in vocab.keys():
      i = vocab[t]
    else:   # if t is [UNK]
      i = 1
    output.append(i)
  return output

さらに、クラスFullTokenizer(object)を以下のように書き換えて、WordpieceTokenizerではなくMecabTokenizerを使うようにします。他の日本語事前学習済みBERTだとMeCabやJuman++で形態素に分割した後でさらにWordpieceTokenizerを適用するものがありますが、ストックマークのモデルではMecabのみを用いています。

bert/tokenization.py

class FullTokenizer(object):
  """Runs end-to-end tokenziation."""

  def __init__(self, vocab_file, do_lower_case=True):
    self.vocab = load_vocab(vocab_file)
    self.inv_vocab = {v: k for k, v in self.vocab.items()}
    #self.basic_tokenizer = BasicTokenizer(do_lower_case=do_lower_case)
    #self.wordpiece_tokenizer = WordpieceTokenizer(vocab=self.vocab)
    # use Mecab
    self.mecab_tokenizer = MecabTokenizer()

  def tokenize(self, text):
    split_tokens = []
    # for token in self.basic_tokenizer.tokenize(text):
    # use Mecab
    for token in self.mecab_tokenizer.tokenize(text):
    	split_tokens.append(token)

    return split_tokens

  def convert_tokens_to_ids(self, tokens):
    #return convert_by_vocab(self.vocab, tokens)
    return convert_tokens_to_ids(self.vocab, tokens)

  def convert_ids_to_tokens(self, ids):
    return convert_by_vocab(self.inv_vocab, ids)

最後にBasicTokenizerを継承したクラスMecabTokenizerを追記します。

bert/tokenization.py

class MecabTokenizer(BasicTokenizer):
  def __init__(self):
    import MeCab
    path = "-d /usr/lib/x86_64-linux-gnu/mecab/dic/mecab-ipadic-neologd"
    self._mecab = MeCab.Tagger(path)

  def tokenize(self, text):
    """Tokenizes a piece of text."""
    text = convert_to_unicode(text.replace(' ', ''))
    text = self._clean_text(text)

    mecab_result = self._mecab.parseToNode(text)
    split_tokens = []
    while mecab_result:
    	if mecab_result.feature.split(",")[0] != 'BOS/EOS':
	        split_tokens.append(mecab_result.surface)
    	mecab_result = mecab_result.next

    output_tokens = whitespace_tokenize(" ".join(split_tokens))
    print(split_tokens)
    return output_tokens

これで日本語事前学習済みのBERTを使う準備ができましたが、extract_features.pyは入力した文章のすべてのトークンに対する埋め込みベクトルを格納したファイルを出力しますので、文章数が多いと出力ファイルの容量が大きくなってしまいます。今回のタスクでは各トークンの埋め込みベクトルではなく、文章の埋め込みベクトルのみを使いますので、それのみを出力するようにextract_features.pyを書き換えることにします。文章の埋め込みベクトルとしては、ELMoと同様にすべてのトークンのベクトルの平均を使います。また、埋め込みベクトルを取り出す層についても、元のコードでは指定した層のベクトルをそれぞれ出力するのですが、指定した層すべてのベクトルの平均を取るように変更します。出力はnumpyのnpy形式で保存するので、ヘッダー部分に

bert/extract_features.py

import numpy as np

を書き加えます。さらに、main関数の最後の部分、input_fn = input_fn_builder(の行以降を以下のように変更します。

bert/extract_features.py

  input_fn = input_fn_builder(
      features=features, seq_length=FLAGS.max_seq_length)

  arr = []
  for result in estimator.predict(input_fn, yield_single_examples=True):
    cnt = 0
    for (i, token) in enumerate(feature.tokens):
      for (j, layer_index) in enumerate(layer_indexes):
        layer_output = result["layer_output_%d" % j]
        if token != '[CLS]' and token != '[SEP]':
          if cnt == 0:
            averaged_emb = np.array(layer_output[i:(i + 1)].flat)
          else:
            averaged_emb += np.array(layer_output[i:(i + 1)].flat)
          cnt += 1
    averaged_emb /= cnt
    arr += [averaged_emb]
  np.save(FLAGS.output_file, arr)

これで準備が整いましたので、ノートブックELMo_BERT.ipynbで以下のコマンドを実行します。埋め込みベクトルを取り出す層としては最終層を除いた後半の５層を使うこととします。

ELMo_BERT_embedding.ipynb

# BERT実行 dev
!python ./bert/extract_features.py \
  --input_file=dev_text.txt \
  --output_file=dev_bert_embeddings.npy \
  --vocab_file=./BERT_base_stockmark/vocab.txt \
  --bert_config_file=./BERT_base_stockmark/bert_config.json \
  --init_checkpoint=./BERT_base_stockmark/output_model.ckpt \
  --layers -2,-3,-4,-5,-6

ELMo_BERT_embedding.ipynb

# BERT実行 test
!python ./bert/extract_features.py \
  --input_file=test_text.txt \
  --output_file=test_bert_embeddings.npy \
  --vocab_file=./BERT_base_stockmark/vocab.txt \
  --bert_config_file=./BERT_base_stockmark/bert_config.json \
  --init_checkpoint=./BERT_base_stockmark/output_model.ckpt \
  --layers -2,-3,-4,-5,-6

USE

USEは日本語を含めた多言語で学習したモデルの軽量版(Multilingual, CNN版, v3)を用います。使用方法は前回の記事で説明した通りです。TensorFlowはバージョン2.x系を使用します。

USE_embedding.ipynb

import tensorflow_hub as hub
import tensorflow_text
import tensorflow as tf

from google.colab import drive
drive.mount('/content/drive')

USE_embedding.ipynb

%cd /content/drive/'My Drive'/anomaly_detection

use_url = 'https://tfhub.dev/google/universal-sentence-encoder-multilingual/3'
embed = hub.load(use_url)

def get_use_embeddings(texts, batch_size=100):
    length = len(texts)
    n_loop = int(length / batch_size) + 1
    embeddings = use_embedder(texts[: batch_size])
    for i in range(1, n_loop):
        arr = use_embedder(texts[batch_size*i: min(batch_size*(i+1), length)])
        embeddings = tf.concat([embeddings, arr], axis=0)
    return np.array(embeddings)

USE_embedding.ipynb

dev_use_embeddings = get_use_embeddings(dev_arr[:, 1])
test_use_embeddings = get_use_embeddings(test_arr[:, 1])
print(dev_use_embeddings.shape, test_use_embeddings.shape)
# (4159, 512) (2060, 512)

np.save('dev_use_embeddings.npy', dev_use_embeddings)
np.save('test_use_embeddings.npy', test_use_embeddings)

異常検知モデル開発

これで３つのモデルの埋め込みベクトルが用意できましたので、異常検知モデルを構築して精度の比較を行います。
ノートブックanomaly_detection.ipynb上で、開発・テストデータと埋め込みベクトルデータを読み込みます。

anomaly_detection.ipynb

import numpy as np
from scipy.stats import chi2

from google.colab import drive
drive.mount('/content/drive')

anomaly_detection.ipynb

%cd /content/drive/'My Drive'/anomaly_detection
dev_arr = np.load('dev.npy')
test_arr = np.load('test.npy')
dev_elmo_embeddings = np.load('dev_elmo_embeddings.npy')
test_elmo_embeddings = np.load('test_elmo_embeddings.npy')
dev_bert_embeddings = np.load('dev_bert_embeddings.npy')
test_bert_embeddings = np.load('test_bert_embeddings.npy')
dev_use_embeddings = np.load('dev_use_embeddings.npy')
test_use_embeddings = np.load('test_use_embeddings.npy')

print(dev_elmo_embeddings.shape, test_elmo_embeddings.shape)
# (4159, 1024) (2060, 1024)
print(dev_bert_embeddings.shape, test_bert_embeddings.shape)
# (4159, 768) (2060, 768)
print(dev_use_embeddings.shape, test_use_embeddings.shape)
# (4159, 512) (2060, 512)

異常検知モデルの詳細については前回の記事を参照してください。モデルのコードを以下のクラスDirectionalAnomalyDetectionにまとめました。

anomaly_detection.ipynb

class DirectionalAnomalyDetection:
  def __init__(self, dev_embeddings, test_embeddings, dev_arr):
    self.dev_embeddings = self.normalize_arr(dev_embeddings)
    self.test_embeddings = self.normalize_arr(test_embeddings)
    self.dev_arr = dev_arr
    self.mu, self.anom  = self.calc_anomality(self.dev_embeddings)
    self.mhat, self.shat = self.calc_chi2params(self.anom)
    print("mhat: {:.3f}".format(self.mhat))
    print("shat: {:.3e}".format(self.shat))
    self.anom_test = None

  def normalize_arr(self, arr):  
    norm = np.apply_along_axis(lambda x: np.linalg.norm(x), 1, arr) # norm of each vector
    normed_arr = arr / np.reshape(norm, (-1,1))
    return normed_arr

  def calc_anomality(self, embeddings):
    mu = np.mean(embeddings, axis=0)
    mu /= np.linalg.norm(mu)
    anom = 1 - np.inner(mu, embeddings)
    return mu, anom

  def calc_chi2params(self, anom):
    anom_mean = np.mean(anom)
    anom_mse = np.mean(anom**2) - anom_mean**2
    mhat = 2 * anom_mean**2 / anom_mse
    shat = 0.5 * anom_mse / anom_mean
    return mhat, shat

  def decide_ath_by_alpha(self, alpha, x_ini, max_ite=100, eps=1.e-12):
    # Newton's method
    x = x_ini
    for i in range(max_ite):
      xnew = x - (chi2.cdf(x, self.mhat, loc=0, scale=self.shat) 
            - (1 - alpha)) / chi2.pdf(x, self.mhat, loc=0, scale=self.shat)
      if abs(xnew - x) < eps:
        print("iteration: ", i+1)
        break
      x = xnew
    print("ath: {:.4f}".format(x))
    return x

  def decide_ath_by_labels(self, x_ini, max_ite=100, eps=1.e-12):
    anom0 = self.anom[self.dev_arr[:, 0] == '0.0']
    anom1 = self.anom[self.dev_arr[:, 0] == '1.0']
    mhat0, shat0 = self.calc_chi2params(anom0)
    mhat1, shat1 = self.calc_chi2params(anom1)
    ath = self._find_crossing_point(mhat0, shat0, mhat1, shat1, x_ini, max_ite, eps)
    print("ath: {:.4f}".format(ath))
    return ath

  def _find_crossing_point(self, mhat0, shat0, mhat1, shat1, x_ini, max_ite, eps):
    # Newton's method
    x = x_ini
    for i in range(max_ite):
      xnew = x - self._crossing_func(x, mhat0, shat0, mhat1, shat1)
      if abs(xnew - x) < eps:
        print("iteration: ", i+1)
        break
      x = xnew
    return x

  def _crossing_func(self, x, mhat0, shat0, mhat1, shat1):
    chi2_0 = chi2.pdf(x, mhat0, loc=0, scale=shat0)
    chi2_1 = chi2.pdf(x, mhat1, loc=0, scale=shat1)
    nume = x * (chi2_0 - chi2_1)
    deno = (mhat0 * 0.5 - 1 - x / shat0 * 0.5) * chi2_0 \
          -(mhat1 * 0.5 - 1 - x / shat1 * 0.5) * chi2_1
    return nume / deno

  def predict(self, ath):
    self.anom_test = 1 - np.inner(self.mu, self.test_embeddings)
    predict_arr = (self.anom_test > ath).astype(np.int)
    return predict_arr

３種類のデータそれぞれに対してモデルインスタンスを作成します。異常度の分布をカイ２乗分布でフィッティングした時のパラメータ$\hat{m}$と$\hat{s}$が出力されます。

anomaly_detection.ipynb

#-- ELMo
dad_elmo = DirectionalAnomalyDetection(dev_elmo_embeddings, test_elmo_embeddings, dev_arr)
# mhat: 6.813
# shat: 3.011e-03

#-- BERT
dad_bert = DirectionalAnomalyDetection(dev_bert_embeddings, test_bert_embeddings, dev_arr)
# mhat: 20.358
# shat: 6.912e-03

#-- USE
dad_use = DirectionalAnomalyDetection(dev_use_embeddings, test_use_embeddings, dev_arr)
# mhat: 36.410
# shat: 1.616e-02

どのモデルに対しても有効次元$\hat{m}$は実際の埋め込みベクトルの次元よりかなり小さくなっています。実際の埋め込みベクトルの次元（ELMo: 1024, BERT: 768, USE: 512）が大きいほど逆に有効次元が小さくなっているのが興味深いです。
こうしてパラメータを決定したカイ２乗分布を実際の異常度の分布と共にプロットすると以下のようになります。異常次元の小ささを反映して、BERTと特にELMoでは異常度の分布のピークの値も小さくなっています。また、カイ2乗分布の当てはまりもELMoとBERTでは比較的良くないのが見て取れます。

次に異常度に対する閾値$a_\text{th}$を求めます。前回の記事では、前もって決めた誤報率$\alpha$に対して式
$$1-\alpha = \int_0^{a_\text{th}} \! dx, \chi^2 (x|\hat{m}, \hat{s})$$を解くことで閾値$a_\text{th}$を決定しました。この方法は上で定義したクラスのメソッドdecide_ath_by_alphaに実装されています。結果だけを示しますが、この方法では３つのどのモデルでもテストデータに対する分類精度はかなり低くなってしまいました。誤報率の値は前回と同じく0.01としています。

	Accuracy	Precision	Recall	F1
ELMo	0.568	0.973	0.141	0.246
BERT	0.580	0.946	0.170	0.288
USE	0.577	0.994	0.155	0.269

なぜこんなに精度が悪いのかを見るために、テストデータに対する異常値のヒストグラムを正常データと異常データに分けてプロットします。左の山が正常データ、右の山が異常データです。赤い縦線がdecide_ath_by_alphaで求めた閾値$a_\text{th}$を表しています。ここではUSEに対する結果だけを示していますが、他のモデルの結果も同様です。正常データと異常データはある程度は分離できていますが、閾値が大きすぎるためにほとんどのデータを正常データとして分類してしまっていることがわかります。

前回の記事で示した結果では、異常度の分布に対するカイ２乗分布の当てはまりが良く、正常データと異常データの分布がかなりはっきりと分かれていたので、上述の方法で決めた閾値で高い精度が得られていました。しかし、今回のデータに対しては、正常データと異常データ分布の重なりが比較的大きいために、うまくいきませんでした。

さて、別の方法で閾値$a_\text{th}$を決定したいのですが、教師なし学習の設定ではこれ以上の手立てはなさそうです。そこで、開発データの正常/異常フラグが既知であるという教師あり学習の設定で閾値を決めることにします。
具体的には以下の手順をとります。
１．開発データの異常度を正常データに対するものと異常データに対するものに分ける。
２．２つに分けた異常度の分布それぞれをカイ２乗分布でフィッティングする。
３．２つのカイ２乗分布の交点を異常度の閾値とする。
USEで計算した異常度に対してこの手順を踏んだ結果を下図に示します。黒と青の実線が２つのカイ２乗分布を、赤い点線が閾値を表しています。

閾値を求めるこの手順はクラスDirectionalAnomalyDetectionのメソッドdecide_ath_by_labelsに実装されていますので、以下のように実行します。

anomaly_detection.ipynb

#-- ELMo
ath_elmo = dad_elmo.decide_ath_by_labels(x_ini=0.02)
# iteration:  5
# ath: 0.0312

#-- BERT
ath_bert = dad_bert.decide_ath_by_labels(x_ini=0.2)
# iteration:  6
# ath: 0.1740

#-- USE
ath_use = dad_use.decide_ath_by_labels(x_ini=0.8)
# iteration:  5
# ath: 0.7924

求めた閾値を使ってテストデータに対する予測をして精度を評価します。

anomaly_detection.ipynb

# 予測
predict_elmo = dad_elmo.predict(ath_elmo)
predict_bert = dad_bert.predict(ath_bert)
predict_use = dad_use.predict(ath_use)

# 正解データ
answer = test_arr[:, 0].astype(np.float)

anomaly_detection.ipynb

# 精度評価のための関数
def calc_precision(answer, predict, title, export_fig=False):
  acc = accuracy_score(answer, predict)
  precision = precision_score(answer, predict)
  recall = recall_score(answer, predict)
  f1 = f1_score(answer, predict)
  cm = confusion_matrix(answer, predict)
  print("Accuracy: {:.3f}, Precision: {:.3f}, Recall: {:.3f}, \
        F1: {:.3f}".format(acc, precision, recall, f1))

  plt.rcParams["font.size"] = 18
  cmd = ConfusionMatrixDisplay(cm, display_labels=[0,1])
  cmd.plot(cmap=plt.cm.Blues, values_format='d')
  plt.title(title)
  if export_fig:
      plt.savefig("./cm_" + title + ".png", bbox_inches='tight')
  plt.show()
  return [acc, precision, recall, f1]

anomaly_detection.ipynb

# 精度評価
_ = calc_precision(answer, predict_elmo, title='ELMo', export_fig=True)
_ = calc_precision(answer, predict_bert2, title='BERT', export_fig=True)
_ = calc_precision(answer, predict_use, title='USE', export_fig=True)

結果は以下の通りです。

	Accuracy	Precision	Recall	F1
ELMo	0.916	0.910	0.923	0.917
BERT	0.851	0.847	0.858	0.852
USE	0.946	0.931	0.963	0.947

混同行列

教師あり学習の設定で決めた閾値を用いることで、どのモデルの精度も大きく向上しています。３つのモデルの比較では、USEの精度がどの指標でも一番高く、逆にBERTの精度は一番低くなっています。

最後に、各モデルに対してテストデータの異常度を正常データと異常データに分けてプロットします。

赤い縦の点線が正解ラベルを使って求めた閾値を表しています。異常データが１％しか含まれない開発データから決定した閾値ですが、確かに正常データと異常データの分布の交わるあたりに位置していることがわかります。BERTに関しては、正常データと異常データの分布の重なりが大きく、それゆえに精度が一番低くなったということが見て取れます。

おわりに

前回の記事のような教師なし学習の設定では異常度の閾値を適切に定めることができなかったため、教師ありの異常検知モデルを作成しました。教師ありといっても正解ラベルは異常度の閾値を決めるのに使っただけなので、精度を左右するのはエンコーダーモデルの良し悪し、すなわち正常データと異常データの異常度の分布をどれだけきれいに分離出来るかに掛かっていました。結果として、一番精度が高かったエンコーダーモデルはUSEでした。コサイン類似度に関するタスクはやはりUSEが得意とするところだということでしょうか。