Universal Sentence Encoderを使って文章の異常検知をする

## 概要
講談社MLPの「異常検知と変化検知」を読んで、何か具体的な問題で試してみたいと思ったので、「方向データの異常検知」を文章の埋め込みベクトルに適用して、文章群に混じった異質な文章を検知できるか試してみました。具体的には、夏目漱石の小説から取った文章群の中に企業の有価証券報告書から取った文章を少数だけ混ぜて、異質なデータである有価証券報告書の文章を検知する機械学習モデルを作成しました。埋め込みベクトル（分散表現）の計算には[Multilingual Universal Sentence Encoder](https://arxiv.org/abs/1907.04307)を用いています。

## 方向データの異常検知
「異常検知と変化検知」（著：井出剛、杉山将）のChapter 7「方向データの異常検知」から必要な事項をまとめます。

正解ラベルの付いていないデータ $\mathcal{D}$ を用いて異常検知モデルを作成するときの基本的な考え方は、データに含まれる異常標本の数が非常に少ないと仮定して、正常データの分布 $p(x|\mathcal{D})$ を推定し、分布から外れている標本を異常とみなすというものです。具体的には異常度を$$a(x^\prime) = - \ln p(x^\prime|\mathcal{D}) \tag{1}$$で定義して、異常度が閾値を超えた時に$x^\prime$を異常と判定します。

データ $\mathcal{D}=\\{\boldsymbol{x}^{(n)}|n=1,\cdots,N\\}$ が単位ベクトル$\boldsymbol{x}^{(n)}$からなる**方向データ**である場合は、データ分布としてフォンミーゼス・フィッシャー分布 ([von Mises-Fisher distribution](https://en.wikipedia.org/wiki/Von_Mises_distribution)) $$\mathcal{M} (\boldsymbol{x} |\boldsymbol{\mu}, \kappa) = \frac{\kappa^{M/2-1}}{(2\pi)^{M/2} I_{M/2-1}(\kappa)} \exp(\kappa \boldsymbol{\mu}^\top \boldsymbol{x}) \tag{2}$$を考えるのが自然です。分布は2つのパラメータ、平均方向$\boldsymbol{\mu}$と集中度$\kappa$を持ちます。平均方向$\boldsymbol{\mu}$も単位ベクトルです。ここで$M$はベクトル$\boldsymbol{x}$の次元、$I_{M/2-1}(\kappa)$は第1種変形ベッセル関数です。分布の係数部分は正規化条件$$\int_{S_M}\\! d\boldsymbol{x}\, \mathcal{M} (\boldsymbol{x} |\boldsymbol{\mu}, \kappa) = 1$$から決まっています。ここで、$S_M$は積分範囲が$M$次元空間の単位球表面であることを表しています。

異常検知モデル作成の流れは非常に単純で、大まかには以下の通りです。
１．与えられたデータ $\mathcal{D}$ を用いて最尤推定などによりパラメータ$\boldsymbol{\mu}$と$\kappa$を決定する。
２．$\mathcal{D}$ に含まれない新たなデータ$\boldsymbol{x}^\prime$に対して異常度$a(\boldsymbol{x}^\prime)$を計算し、閾値を超えるかどうかで異常かどうかを判定。

まず、$\boldsymbol{\mu}$を決定する最尤推定をします。条件$\boldsymbol{\mu}^\top \boldsymbol{\mu}=1$をラグランジュ係数$\lambda$を用いて取り込むと$$0=\frac{\partial}{\partial \boldsymbol{\mu}} \left\\{ \ln \prod_{n=1}^N \mathcal{M} (\boldsymbol{x}^{(n)} |\boldsymbol{\mu}, \kappa) -\lambda \boldsymbol{\mu}^\top \boldsymbol{\mu} \right\\}=\kappa \sum_{n=1}^N \boldsymbol{x}^{(n)} -2\lambda \boldsymbol{\mu}$$が得られます。条件$\boldsymbol{\mu}^\top \boldsymbol{\mu}=1$を使って係数$\lambda$を消去すれば、$\boldsymbol{\mu}$をデータ $\mathcal{D}$ から決定する式$$\hat{\boldsymbol{\mu}}=\frac{\boldsymbol{m}}{\sqrt{\boldsymbol{m}^\top \boldsymbol{m}}} \hspace{10pt} \text{where} \hspace{10pt} \boldsymbol{m}=\frac{1}{N}\sum_{n=1}^N \boldsymbol{x}^{(n)} \tag{3}$$が得られます。これはデータベクトルの平均方向が分布の中心に一致するという結果で直感的にも納得のいくものでしょう。

次に$\kappa$を決定したいのですが、分布(2)はベッセル関数を通じて複雑に$\kappa$に依存しているので、最尤推定をしても方程式を解析的に解くことはできません。そこで、モーメント法を用いて$\kappa$を推定します。ここでいうモーメント法とは、仮定した分布から何らかの量のモーメント（期待値や2乗期待値など）を計算し、それらをデータ標本から計算したモーメントと等値することで、仮定した分布のパラメータを決定する方法です。ここでは、異常度に対してモーメント法を適用します。まず、フォンミーゼス・フィッシャー分布(2)を異常度(1)に代入すると $a(\boldsymbol{x}^\prime) = \text{const} -\kappa \hat{\boldsymbol{\mu}}^\top \boldsymbol{x}^\prime$ が得られます。constは$\boldsymbol{x}^\prime$に依存しない定数です。さらに、異常度に定数をかけたり足したりする操作は閾値の再定義で吸収できる事実を使って、$$a(\boldsymbol{x}^\prime) = 1 -\hat{\boldsymbol{\mu}}^\top \boldsymbol{x}^\prime \tag{4}$$とすることができます。単位ベクトルどうしの内積の取りうる値の範囲は$[-1,1]$ですので、異常度の取りうる値の範囲は$0\leq a\leq 2$です。$\boldsymbol{x}$がフォンミーゼス・フィッシャー分布(2)に従うとき、異常度$a$の確率分布は$$p(a)=\int_{S_M}\\! d\boldsymbol{x}\, \delta\left(a-(1 -\hat{\boldsymbol{\mu}}^\top \boldsymbol{x})\right) \mathcal{M} (\boldsymbol{x} |\hat{\boldsymbol{\mu}}, \kappa) $$と書き下せます。デルタ関数を用いて積分を実行し、データ点の平均方向$\hat{\boldsymbol{\mu}}$周りのばらつきは小さい、すなわち$a\ll 1$であると仮定すると、近似的に$$p(a) \simeq c\, a^{(M-1)/2-1}\exp(-\kappa a)$$が得られます。$c$は$a$に依らない定数です。この分布はカイ2乗分布$$\chi^2 (x|m,s)=\frac{1}{2s\Gamma(m/2)}\left(\frac{x}{2s}\right)^{m/2-1}\exp\left(-\frac{x}{2s}\right) \tag{5}$$で$m=M-1$, $s=1/(2\kappa)$としたものに一致します。異常度$a$の従う分布がわかりましたので、1次と2次のモーメントを以下のように計算することができます[^1]。$$\langle a\rangle=\int_0^\infty \\! da\, a\, \chi^2(a|m,s)=ms$$ $$\langle a^2\rangle=\int_0^\infty \\! da\, a^2\, \chi^2(a|m,s)=m(m+2)s^2 .$$一方で、モーメントはデータ標本$\boldsymbol{x}^{(n)}$から$$\langle a\rangle \simeq \frac{1}{N}\sum_{n=1}^N a(\boldsymbol{x}^{(n)}) \ , \hspace{10pt} \langle a^2\rangle \simeq \frac{1}{N}\sum_{n=1}^N a(\boldsymbol{x}^{(n)})^2$$と評価することもできます。2つの表現を等値することで、パラメータの推定値$$\hat{m}=\frac{2\langle a\rangle^2}{\langle a^2\rangle -\langle a\rangle^2}  \ , \hspace{10pt} \hat{s} =\frac{\langle a^2\rangle -\langle a\rangle^2}{2\langle a\rangle} \tag{6}$$が得られます。こうして得られた$\hat{m}$は理想的には$M-1$（$M$はベクトルの次元）に等しいはずですが、後で示すように実際は$M$よりかなり小さくなります。これはデータの主要なばらつきが一部の変数に集中しているからで、その意味で$\hat{m}$は有効次元と呼ばれます。はじめの目的だった$\kappa$の推定値は関係式$\hat{\kappa}=1/(2\hat{s})$から求められますが、異常度を(4)式のように$\kappa$を含まない形で再定義しているので、実は$\hat{\kappa}$の値を直接使う必要はもはやありません。異常度の閾値を決めるプロセスで$\hat{m}$と$\hat{s}$を用いるので、そこで間接的に$\hat{\kappa}$の値を使ったことになります。

さて、異常度$a$の従う分布が$\chi^2 (a|\hat{m}, \hat{s})$と推定できましたので、新たなデータ$\boldsymbol{x}^\prime$に対して異常度を式(4)を用いて計算し、その値が分布から外れていれば異常データであると判定できます。判定基準としては、誤報率$\alpha$を何らかの経験則により決めて、$$1-\alpha = \int_0^{a_\text{th}} \\! dx\, \chi^2 (x|\hat{m}, \hat{s}) \tag{7}$$を解くことで、異常度の閾値$a_\text{th}$を求めることができます[^2]。上の式から明らかなように、誤報率$\alpha$は確率分布の端の何パーセントに入ったら異常データという基準を与えるものです。以下で示す実装では、ニュートン法を用いてこの式を解いています。

[^1]: 異常度$a$の取りうる値の範囲は$0\leq a\leq 2$である上に、カイ2乗分布は$a\ll 1$という仮定の下で導かれたものですが、被積分関数は大きな$a$に対しては非常に小さな値を取るので、問題なく積分範囲を$+\infty$まで拡張することができます。
[^2]: 式(7)の右辺はガンマ関数と不完全ガンマ関数を用いて$\gamma(a_\text{th}/2, \hat{m}/2)/\Gamma(\hat{m}/2)$と書くこともできます。

## Universal Sentence Encoder
Universal Sentence Encoder (USE)は任意の長さの文章を512次元のベクトルに変換するエンコーダーモデルです。自然言語処理の様々なタスクで事前学習されています。BERTなどの汎用言語モデルと違って、正解ラベル付きの教師あり学習により事前学習されているので、事前学習と同種のタスクに対してはファインチューニングなしでも高い性能を発揮します。特に、文章の類似度を判定するタスクはコサイン類似度をターゲットに学習されているので、コサイン類似度を使ったタスクと相性が良いです。今回の異常検知モデルはベクトルの内積から異常度を計算しますので、USEに向いたタスクであると期待できます。

2018年にGoogleが発表した[最初のバージョン](https://arxiv.org/abs/1803.11175)は英語の文章にのみ対応していましたが、2019年に発表された[Multilingual USE](https://arxiv.org/abs/1907.04307)は日本語を含む16言語に対応しています。モデルが16種類あるわけではなく、単一のモデルが16言語に対応しています。多言語バージョンは同一言語内での学習タスクと異言語間での翻訳タスクで事前学習されています。英語以外の学習データを水増しするために、なんと英語のデータセットをGoogle翻訳したそうです。

モデルのアーキテクチャは（論文にもあまり詳細が書かれていないのですが）英語バージョンは[Transformer](https://arxiv.org/abs/1706.03762)を用いたものと、[Deep Averaging Network](https://people.cs.umass.edu/~miyyer/pubs/2015_acl_dan.pdf) (DAN)を用いたものがあります。DANの方は精度は落ちるものの軽量なモデルとなっています。多言語バージョンもアーキテクチャは2種類あるのですが、軽量モデルはDANではなくCNNとなっています。TransformerでもCNNでもトークンレベルの埋め込みベクトルを平均したものを文章の埋め込みベクトルとしているようです。文章のトークン化には[SentencePiece](https://qiita.com/taku910/items/7e52f1e58d0ea6e7859c)が使われており、日本語を入力する場合にも事前に分かち書きをする必要はありません。

いずれのバージョンの学習済みモデルもTensorFlow Hubから使用でき、多言語・Transformer版は[USE-multilingual-large](https://tfhub.dev/google/universal-sentence-encoder-multilingual-large/3)、多言語・CNN版は[USE-multilingual](https://tfhub.dev/google/universal-sentence-encoder-multilingual/3)からアクセスできます。今回の実験では軽量なCNN版を使用します[^3]。また、USEモジュールの出力は初めから単位ベクトルに正規化されています。

[^3]: 実はTransformer版も使ってみたのですが、CNN版よりもわずかに精度が悪くなってしまいました。


## 実装
### 環境
Ubuntu 18.04.4 LTS (docker on windows)
Python 3.7.6
tensorflow 2.1.0
tensorflow-hub 0.8.0
tensorflow-text 2.1.1

(2021/6/21追記) 
Githubにコードを上げました。（各ライブラリのバージョンは上記とは異なっています。）
https://github.com/jovyan70/sentence-anomaly-detection

### データ準備
夏目漱石の小説の文章は後期3部作「彼岸過迄」「行人」「こころ」を使用することとします。[青空文庫](http://www.aozora.gr.jp/)からBeautifulSoupを用いてスクレイピングします。以下のサイトを参考にしました。

- [青空文庫の自然言語処理をサクッとやる](https://qiita.com/yutopeer/items/92d75865a72c60ae7285)
- [BeautifulSoupの文字化けが止まらない時の解消方法](https://lets-hack.tech/programming/languages/python/beautifulsoup-mojibake/)

```python
import json
import requests
from bs4 import BeautifulSoup
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import random
from sklearn.model_selection import train_test_split

url_list = ['https://www.aozora.gr.jp/cards/000148/files/773_14560.html',  # 彼岸過迄
            'https://www.aozora.gr.jp/cards/000148/files/765_14961.html',  # 行人
            'https://www.aozora.gr.jp/cards/000148/files/775_14942.html',  # こころ
           ]
```
```python
def fetch_sentences(aozora_url):
    res = requests.get(aozora_url)
    soup = BeautifulSoup(res.content, 'html.parser', from_encoding='utf-8')
    main_text = soup.find("div", class_='main_text')

    #ルビが振ってあるのを削除
    for yomigana in main_text.find_all(["rp","h4","rt"]):
        yomigana.decompose()
        
    sentences = [line.strip() for line in main_text.text.strip().splitlines()]
    return sentences

souseki_texts = []
for url in url_list:
    souseki_texts += fetch_sentences(url)

print(len(souseki_texts))
# 6102
```
リスト`souseki_texts`の要素は、いくつかの文をまとめたパラグラフになっていますので、句点ごとに分割することとします。これは、この後ダウンロードするchABSAデータセットの有価証券報告書の文章と形式を合わせるためです。

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

souseki_texts = split_texts(souseki_texts)

print(len(souseki_texts))
# 21011
```
さらに、短すぎる文章と長すぎる文章を除きます。リストには空白も含まれていますが、この操作で同時に除去されます。

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

souseki_texts = filter_by_length(souseki_texts)

print(len(souseki_texts))
# 12221
```
文章の前処理（正規化）は特に行いません。

次に、[chABSA-dataset](https://github.com/chakki-works/chABSA-dataset)をダウンロードして有価証券報告書の文章を採取します。chABSA-datasetは上場企業の有価証券報告書をベースに作成されたデータセットで、各文章の部分部分にネガティブ・ポジティブのラベルがついたものです。詳しい使用例については

- [BERTを用いたネガポジ分類機の作成](https://qiita.com/ysiny/items/b01250228e0c5cc0e647)

が参考になります。今回はラベルは使用せず有価証券報告書の文章のみを使います。

```bash
!wget https://s3-ap-northeast-1.amazonaws.com/dev.tech-sketch.jp/chakki/public/chABSA-dataset.zip && \
 unzip chABSA-dataset.zip
```
```python
data_path = Path("chABSA-dataset")
path_list = [p for p in data_path.glob("*.json")]

chabsa_texts = []
for p in 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
```
漱石の文章と同様に、短すぎる文章と長すぎる文章を除きます。

```python
chabsa_texts = filter_by_length(chabsa_texts)
print(len(chabsa_texts))
# 5148
```
これで2種類の文章群の準備ができましたので、2つを混ぜてデータセットを作ります。まず、漱石にラベル0、有価証券報告書にラベル1を付与します。

```python
souseki_length = len(souseki_texts)
souseki_arr = np.hstack([np.reshape(np.zeros(souseki_length), (-1, 1)), np.reshape(souseki_texts, (-1, 1))])

chabsa_length = len(chabsa_texts)
chabsa_arr = np.hstack([np.reshape(np.ones(chabsa_length), (-1, 1)), np.reshape(chabsa_texts, (-1, 1))])
```
データセットはモデル開発用(dev)とテスト用(test)の2つ用意します。漱石の文章の8割を開発用に、残り2割をテスト用に使います。開発データには漱石の文章の1%の量だけ有価証券報告書の文章を混ぜます。テスト用に関しては、精度評価をしやすいように漱石の文章と有価証券報告書の文章は同数とします。

```python
num_s_dev = int(souseki_length * 0.8)
num_c_dev = int(num_s_dev * 0.01)
num_s_test = souseki_length - num_s_dev
num_c_test = num_s_test

print("開発データ　 漱石:{}, 有報:{}".format(num_s_dev, num_c_dev))
print("テストデータ 漱石:{}, 有報:{}".format(num_s_test, num_c_test))
# 開発データ　 漱石:9776, 有報:97
# テストデータ 漱石:2445, 有報:2445

s_dev, s_test = train_test_split(souseki_arr, train_size=num_s_dev)
c_dev, c_test = train_test_split(chabsa_arr, train_size=num_c_dev)
c_test, _ = train_test_split(c_test, train_size=num_c_test)

dev_arr = np.vstack([s_dev, c_dev])
np.random.shuffle(dev_arr)
test_arr = np.vstack([s_test, c_test])
np.random.shuffle(test_arr)

print(dev_arr.shape, test_arr.shape)
# (9873, 2) (4890, 2)
```
テストデータの最初のいくつかをサンプルしてみると、以下のようになります。
> 1: 可塑剤の出荷数量は増加したものの、原料価格の下落に伴う販売価格の引き下げにより、売上高は46億３百万円となり、前連結会計年度に比べ１億25百万円の減収（同2.7％減）となりました
0: その或ものは水の色を離れない蒼い光を鱗に帯びて、自分の勢で前後左右に作る波を肉の裏に透すように輝やいた
0: 彼はやむを得ず、手紙で用を弁ずる事にした
1: 健康ビジネス部門については、Ｅウール毛布シリーズが順調に売上を伸ばしたものの、その他の健康寝具関係での対前年の反動減が響き、売上高は前期を下回った
1: また製紙薬剤関連原料は前年度実績を下回りましたものの、炭素製品は堅調に推移し、インキ用原料は前年度並みとなりました
0: 未亡人は私の身元やら学校やら専門やらについて色々質問しました
0: 要するに、勝っても負けても、炬燵にあたって、将碁を差したがる男であった
0: お兼さんだけは依然として元の席を動かなかった
1: 当事業年度の経営成績は次のとおりであります
0: 彼はそう云って、すたすた山の方へ歩いて行ったそうです

人間の目で見ると容易に夏目漱石か有価証券報告書かを判別できそうですね。それでは、異常検知モデルを構築して、モデルが2種類の文章を判別をできるかどうか見ていきましょう。

### 異常検知モデル開発
まず初めに、Universal Sentence Encoder (USE)を用いて、文章を埋め込みベクトルに変換します。必要なライブラリをインポートし、TensorFlow Hubからモデル（[Multilingual, CNN版](https://tfhub.dev/google/universal-sentence-encoder-multilingual/3)）をロードします。

```python
import tensorflow_hub as hub
import tensorflow_text
import tensorflow as tf

use_url = 'https://tfhub.dev/google/universal-sentence-encoder-multilingual/3'
embed = hub.load(use_url)
```
USEモデルの使い方は単純で、ロードした関数`embed`に文章のリストを与えれば、tf.Tensor形式の埋め込みベクトル（単位ベクトル）のリストを返してくれます。全ての文章を一括で渡すとメモリの使用量が大きくなるので、バッチを分割して計算します。

```python
def get_embeddings(texts, batch_size=100):
    length = len(texts)
    n_loop = int(length / batch_size) + 1
    embeddings = embed(texts[: batch_size])
    for i in range(1, n_loop):
        arr = embed(texts[batch_size*i: min(batch_size*(i+1), length)])
        embeddings = tf.concat([embeddings, arr], axis=0)
    return np.array(embeddings)
```
```python
dev_embeddings = get_embeddings(dev_arr[:, 1])
test_embeddings = get_embeddings(test_arr[:, 1])
print(dev_embeddings.shape, test_embeddings.shape)
# (9873, 512) (4890, 512)
```
こうして得られた、`dev_embeddings`を式(3)に代入して、平均方向の推定値$\hat{\boldsymbol{\mu}}$を求めます。

```python
mu = np.mean(dev_embeddings, axis=0)
mu /= np.linalg.norm(mu)
print(mu.shape)
# (512,)
```
次に、式(4)により開発データの各標本に対して異常度を計算し、それを式(6)に代入することで、異常度の従うカイ2乗分布のパラメータ$\hat{m}$と$\hat{s}$を推定します。

```python
anom = 1 - np.inner(mu, dev_embeddings)
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
print("mhat: {:.1f}".format(mhat))
print("shat: {:.3e}".format(shat))
# mhat: 111.7
# shat: 5.406e-03
```
すでに述べたように、$\hat{m}$は理想的には$M-1$に一致するはずですが、実際の値は111.7で$M-1=511$よりかなり小さくなっています。

求めたパラメータに対して、カイ2乗分布の確率密度関数PDF(5)と累積分布関数CDF（式(7)の右辺）をプロットすると以下のようになります。左図には実際の異常度の分布もヒストグラムでプロットしています。確かにカイ2乗分布でフィット出来ていることが見て取れます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613488/6b542c67-868c-4258-b99a-6a0af11a45b1.png", height=250><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613488/5e754e82-88b1-f8a6-d819-ff98589cc45c.png", height=250>



最後に式(7)を解いて閾値$a_\text{th}$を求めれば異常検知モデルが完成します。そのためには誤報率$\alpha$を設定する必要があります。誤報率は何らかのヒューリスティクスによって決める必要がありますが、今回は開発データには異常データが約1%含まれることを知っていますので、$\alpha=0.01$と設定します。式(7)はニュートン法で以下のように解くことができます。

```python
from scipy.stats import chi2

alpha = 0.01
x = 0.8  # initial value
for i in range(100):
    xnew = x - (chi2.cdf(x, mhat, loc=0, scale=shat) - (1 - alpha)) / chi2.pdf(x, mhat, loc=0, scale=shat)
    if abs(xnew - x) < 1.e-12:
        print("iteration: ", i+1)
        break
    x = xnew
ath = x
print("ath: {:.4f}".format(ath))
# iteration:  5
# ath: 0.8077
```
### 精度評価
モデルが完成しましたので、テストデータを使って精度評価を行います。
まず、開発データから評価した平均方向$\hat{\boldsymbol{\mu}}$を使って、テストデータに対する異常度を計算します。

```python
anom_test = 1 - np.inner(mu, test_embeddings)
```
得られた異常度を閾値と比較して、閾値より大きければ異常標本として判定します。

```python
predict = (anom_test > ath).astype(np.int)
```
予測値を正解データと比較することで各種精度指標を計算します。混同行列の表示には[ConfusionMatrixDisplay](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.ConfusionMatrixDisplay.html)を利用しました。

```python
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, confusion_matrix, ConfusionMatrixDisplay

answer = test_arr[:, 0].astype(np.float)

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.show()
```
結果は以下の通りです。

|Accuracy|Precision|Recall|F1|
|--:|--:|--:|--:|
|0.990|0.993|0.987|0.990|
混同行列
![cm.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613488/e46815c8-330e-6fa0-9dcf-cb0e845ef67a.png)
単純なモデルにも関わらず思いのほか高い精度が得られました。Accuracy, Precision, Recallがどれも同じくらいの値となっていますので、誤報率$\alpha$の設定値も適切であったと結論できます。

最後に、テストデータに対する異常値のヒストグラムを正常データと異常データに分けてプロットします。左の山が正常データ、右の山が異常データです。黒い実線は開発データから推定した$\hat{m}$と$\hat{s}$に対するカイ2乗分布です。赤の縦線は閾値$a_\text{th}$を表していますので、閾値のあたりを境に2つの分布が良く分けられていることがわかります。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/613488/168af4ea-6c87-fb99-5e46-39b5cc8b9aca.png", height=300>

## おわりに
ちょっとした思い付きで計算して記事を書いてみたのですが、個人的には思いのほか面白くなって記事も長くなってしまいました。今回試したのは人間が見れば簡単に見分けられるような簡単なタスクでしたが、正直ここまで精度が高くなるとは思っておらず驚きました。異常検知モデル自体は単純なものですので、USEの性能の高さが精度につながったのではないかと思います。その意味では、他のエンコーダーモデルでも同じ計算をして比較することも必要だったと感じています。余裕があれば、後から計算して書き足したいと思います。