Hayami Kaito Web
~/tech

異常検知のしくみ:外れ値・マハラノビス距離から LOF・Isolation Forest・再構成誤差まで

published

工場の機械の故障の予兆、クレジットカードの不正利用、サーバーへの攻撃、製品の傷など、「いつもと違う」データを見つけたい場面はたくさんあります。これが 異常検知です。

異常検知の難しさは、異常の例がほとんど集まらないことです。機械はめったに壊れませんし、不正の手口は次々に変わります。そのため、画像分類 のように「正常」と「異常」の例を大量に集めて分類を学ぶ方法は、そのままでは使えません。かわりに、正常なデータがどんな姿をしているかを学び、そこから外れたものを異常とみなします。

結論を先に書くと、異常検知は次の考え方でできています。

  • 平均から何標準偏差離れているかで外れ値を見つけるのが基本。ただし、異常が多く混じると平均と標準偏差が引っ張られるので、中央値を使う方法もある
  • 1 つずつ見ると普通でも、組み合わせがおかしいデータがある。変数どうしの相関を考えたマハラノビス距離で見つけられる
  • 正常なデータがいくつものまとまりに分かれていると、全体の平均からの距離は役に立たない。近所の点との距離や密度(k 近傍・LOF)、孤立させやすさ(Isolation Forest)で測る
  • 多くのセンサーが連動して動くデータでは、圧縮して戻したときのずれ(再構成誤差)で、連動の崩れを見つけられる
  • 異常はめったにないので、正解率ではなく、適合率と再現率(PR-AUC)で評価する

説明の途中に出てくる計算例は、Python と NumPy で実際に動かした結果です。

環境

計算例を動かした環境です。

項目バージョン
Python3.13
NumPy2.5

異常検知の問題の形

異常検知は、使えるデータによって、次のように分かれます。

使えるデータ方法例
正常と異常の両方に正解がある分類の学習(教師あり)過去の不正利用の記録が大量にある
正常なデータだけがある正常なデータの姿を学び、外れたものを異常とする正常に動いていたときのセンサーの記録
正解のないデータに、異常が少し混じっているデータの大部分から外れたものを異常とする(教師なし)アクセスの記録全体から、怪しいものを探す

異常の例が十分にあれば分類の学習がいちばん強力ですが、ふつうは集まりません。また、過去に見たことのない種類の異常にも気付きたいので、この記事では下の 2 つを扱います。

1 つの値の外れ値

平均から何標準偏差離れているか

いちばん基本的な方法は、値が平均から標準偏差の何倍離れているか(z スコア)を計算し、3 倍を超えたら異常とする方法です。データが正規分布に従うなら、正常な値が平均から 3 標準偏差以上離れる確率は約 0.3% しかありません。

ふだんは 50 度前後で、故障したときだけ 80 度前後になる機械の温度で試します。1,000 回の記録のうち、故障は 10 回です。

import numpy as np

rng = np.random.default_rng(0)
# 1. 1 次元:機械の温度。ふだんは 50 度前後で、故障した 10 回だけ 80 度前後
temp = np.r_[50 + 2 * rng.standard_normal(990), 80 + 3 * rng.standard_normal(10)]
z = (temp - temp.mean()) / temp.std()
print(round(temp.mean(), 2), round(temp.std(), 2), int((np.abs(z) > 3).sum()))

実行結果:

50.2 3.55 10

平均は 50.2 度、標準偏差は 3.55 度で、z スコアが 3 を超えたのは 10 回、つまり故障の 10 回をすべて見つけられました。

異常が多いと平均と標準偏差が引っ張られる

ところが、故障が 1,000 回のうち 100 回に増えると、様子が変わります。上のコードの続きです。

# 異常が多く混じると、平均と標準偏差が引っ張られて見逃す
temp2 = np.r_[50 + 2 * rng.standard_normal(900), 80 + 3 * rng.standard_normal(100)]
z2 = (temp2 - temp2.mean()) / temp2.std()
print(round(temp2.mean(), 2), round(temp2.std(), 2), int((np.abs(z2[900:]) > 3).sum()))
# 中央値と MAD(中央値からの絶対偏差の中央値)を使うと、異常に引っ張られにくい
med = np.median(temp2); mad = 1.4826 * np.median(np.abs(temp2 - med))
z3 = (temp2 - med) / mad
print(round(med, 2), round(mad, 2), int((np.abs(z3[900:]) > 3).sum()), int((np.abs(z3[:900]) > 3).sum()))

実行結果:

52.97 9.13 39
50.4 2.34 100 0

1 行目のとおり、故障の値に引っ張られて、平均は 52.97 度、標準偏差は 9.13 度まで大きくなりました。その結果、80 度前後の故障でも z スコアが 3 を超えにくくなり、100 回のうち 39 回しか見つけられません。異常そのものが、異常を見つけるものさしを狂わせてしまうのです。

2 行目は、平均のかわりに 中央値、標準偏差のかわりに MAD(中央値からの絶対偏差の中央値。1.4826 を掛けると、正規分布の標準偏差と同じ目盛りになる)を使った結果です。中央値は全体の真ん中の値なので、1 割の異常が混じっても 50.4 度とほとんど動かず、ものさしは 2.34 度と、正常な温度のばらつき(2 度)に近いままです。故障 100 回をすべて見つけ、正常な 900 回は 1 回も誤って異常としませんでした。このように、外れ値に引っ張られにくい統計の方法を 頑健な(ロバストな)方法と呼びます。

組み合わせがおかしいデータ

1 つずつ見ると普通の値

複数の値を同時に見るときは、それぞれの値は普通でも、組み合わせがおかしいデータがあります。身長と体重で考えてみます。身長 180 cm も体重 55 kg も珍しくはありませんが、「身長 180 cm で体重 55 kg」は珍しい組み合わせです。上のコードの続きです。

# 2. 2 次元:身長と体重。どちらも普通の値だが、組み合わせがおかしい点
n = 1000
height = 170 + 6 * rng.standard_normal(n)
weight = 65 + 0.9 * (height - 170) + 4 * rng.standard_normal(n)
X = np.c_[height, weight]
odd = np.array([[180.0, 55.0], [160.0, 75.0]])   # 背が高いのに軽い、背が低いのに重い
mu, sd = X.mean(axis=0), X.std(axis=0)
print(np.round((odd - mu) / sd, 2))                # 1 つずつ見ると、どちらも 2 標準偏差以内

実行結果:

[[ 1.74 -1.5 ]
 [-1.65  1.52]]

「180 cm・55 kg」と「160 cm・75 kg」の 2 人は、身長も体重も、z スコアで見ると平均から 2 標準偏差以内に収まっています。1 つずつ見る方法では、異常とはみなされません。

マハラノビス距離:相関を考えた距離

身長と体重には、「背が高い人ほど重い」という相関があります。この相関を考えて距離を測るのが マハラノビス距離です。データの共分散行列の逆行列を使って、データが広がっている向きには距離を短く、広がっていない向きには距離を長く測ります。上のコードの続きです。

cov = np.cov(X.T)
inv = np.linalg.inv(cov)
def mahalanobis(P):
    d = P - mu
    return np.sqrt(np.einsum("ij,jk,ik->i", d, inv, d))
print(round(np.corrcoef(X.T)[0, 1], 2), np.round(mahalanobis(odd), 2), round(float(np.quantile(mahalanobis(X), 0.99)), 2))

実行結果:

0.79 [5.04 4.93] 2.92

身長と体重の相関係数は 0.79 です。2 人のマハラノビス距離は 5.04 と 4.93 で、1,000 人のデータの 99% が収まる距離(2.92)を大きく超えています。データの点の集まりは「右上がりの細長い楕円」の形をしていて、2 人はその楕円の 細い方向に外れているので、遠いと判定されるのです。

マハラノビス距離の 2 乗は、品質管理で ホテリングの T² として使われてきた、多変量の異常検知の基本です。ただし、データ全体が 1 つの楕円(正規分布)の形をしていることを前提にしています。

まとまりがいくつもあるデータ

4 つの方法を比べる

実際のデータは、1 つの楕円の形をしているとは限りません。たとえば、平日と休日でサーバーの負荷のパターンが違えば、正常なデータが 2 つのまとまりに分かれます。

点が密に集まったまとまり(600 点)と、まばらに広がったまとまり(380 点)からなるデータに、異常を 20 点混ぜます。そのうち 10 点は全体から外れた場所に、残りの 10 点は 密なまとまりのすぐ外に置きます。後者は、まばらなまとまりの中の点どうしの間隔と比べれば近い場所にありますが、密なまとまりから見れば明らかに浮いている点です。どれが異常かは教えずに、次の 4 つの方法で、それぞれの点の「異常らしさ」の点数を計算します。

  • マハラノビス距離:全体を 1 つの楕円とみなして、中心からの距離を測る
  • k 近傍距離:10 番目に近い点までの距離。まわりに点が少ないほど大きくなる
  • LOF(局所外れ値因子):自分のまわりの密度が、近所の点のまわりの密度と比べて何分の 1 かを測る。1 なら近所と同じくらいの密度で、大きいほど浮いている
  • Isolation Forest:ランダムに選んだ軸のランダムな位置でデータを切る操作をくり返し、その点が 1 点だけになるまでの回数を数える。外れた点は少ない回数で孤立する。100 本の木で平均する

異常はめったにないので PR-AUC で評価する

異常検知では、異常がとても少ないので、正解率で評価してはいけません。この例では、全部を「正常」と答えるだけで正解率は 98% になってしまいます。そこで、機械学習の評価指標のしくみ で説明した次の指標を使います。

  • ROC-AUC:異常と正常を 1 つずつ選んだとき、異常のほうに高い点数を付けられる確率
  • PR-AUC(平均適合率):点数の高い順に調べていったときの、適合率と再現率の曲線の下の面積。異常が少ないときの性能の差がよく表れる
import numpy as np

rng = np.random.default_rng(0)
# ふだんのデータ:密なまとまり(600 点)と、まばらなまとまり(380 点)
dense = 0.3 * rng.standard_normal((600, 2))
sparse = np.array([6.0, 6.0]) + 1.5 * rng.standard_normal((380, 2))
# 異常 20 点:全体から外れた点 10 個と、密なまとまりのすぐ外にある点 10 個
angle = rng.uniform(0, 2 * np.pi, 10)
far = rng.uniform(-4, 12, (10, 2))
near = np.c_[1.5 * np.cos(angle), 1.5 * np.sin(angle)]
X = np.r_[dense, sparse, far, near]
y = np.r_[np.zeros(980), np.ones(20)]                      # 1 が異常

def roc_auc(score):
    # 異常と正常を 1 つずつ選んだとき、異常のほうが高い点数になる確率
    pos, neg = score[y == 1], score[y == 0]
    return round(float((pos[:, None] > neg[None]).mean()), 3)

def average_precision(score):
    # 点数の高い順に見ていき、異常を見つけるたびの適合率を平均する(PR 曲線の下の面積)
    order = np.argsort(-score)
    hits = y[order]
    precision = np.cumsum(hits) / np.arange(1, len(y) + 1)
    return round(float((precision * hits).sum() / hits.sum()), 3)

def report(name, score):
    top = np.argsort(-score)[:20]
    # 点数の上位 20 点に入った異常の数(全体)と、そのうち密なまとまりのすぐ外の異常の数
    print(f"{name:8s}", roc_auc(score), average_precision(score), int(y[top].sum()), int((top >= 990).sum()))

# 1. マハラノビス距離(全体を 1 つの正規分布とみなす)
mu, inv = X.mean(axis=0), np.linalg.inv(np.cov(X.T))
d = X - mu
report("マハラノビス", np.einsum("ij,jk,ik->i", d, inv, d))

# 2. k 近傍距離:k 番目に近い点までの距離
dist = np.linalg.norm(X[:, None] - X[None], axis=2)
np.fill_diagonal(dist, np.inf)
k = 10
knn = np.sort(dist, axis=1)[:, :k]
report("k近傍距離", knn[:, -1])

# 3. LOF:近所の点の密度と比べて、自分のまわりがどれだけまばらか
neighbors = np.argsort(dist, axis=1)[:, :k]
k_dist = knn[:, -1]
reach = np.maximum(dist[np.arange(len(X))[:, None], neighbors], k_dist[neighbors])
lrd = 1 / reach.mean(axis=1)                               # 局所的な密度
report("LOF", lrd[neighbors].mean(axis=1) / lrd)

# 4. Isolation Forest:ランダムに切っていき、少ない回数で 1 点だけに孤立する点ほど異常
def build(data, depth=0, limit=8):
    if len(data) <= 1 or depth >= limit:
        return len(data)                                   # 葉:残った点の数
    f = rng.integers(2)
    lo, hi = data[:, f].min(), data[:, f].max()
    if lo == hi:
        return len(data)
    cut = rng.uniform(lo, hi)
    side = data[:, f] < cut
    return (f, cut, build(data[side], depth + 1), build(data[~side], depth + 1))

def c(n):
    # n 点が残った葉の先で、さらに切るのに平均で必要な回数の見積もり
    return 2 * (np.log(n - 1) + 0.5772) - 2 * (n - 1) / n if n > 1 else 0

def path_length(x, node, depth=0):
    if not isinstance(node, tuple):
        return depth + c(node)
    f, cut, left, right = node
    return path_length(x, left if x[f] < cut else right, depth + 1)

trees = [build(X[rng.choice(len(X), 256, replace=False)]) for _ in range(100)]
depth = np.array([np.mean([path_length(x, t) for t in trees]) for x in X])
report("iForest", -depth)

実行結果(この環境では 7 秒ほどかかりました):

マハラノビス   0.631 0.303 6 0
k近傍距離    0.956 0.485 8 0
LOF      0.998 0.927 16 9
iForest  0.973 0.578 11 2

左から、ROC-AUC、PR-AUC、点数の上位 20 点に入った異常の数、そのうち密なまとまりのすぐ外の異常の数です。

  • マハラノビス距離は、PR-AUC 0.303 といちばん悪い結果でした。2 つのまとまりを 1 つの大きな楕円とみなすので、楕円の中心付近(2 つのまとまりの間)にある点を正常と判定してしまいます
  • k 近傍距離は ROC-AUC が 0.956 と高くなりましたが、密なまとまりのすぐ外の異常は 1 つも上位 20 点に入りません。そこから 10 番目に近い点までの距離は、まばらなまとまりの中の正常な点と同じくらいだからです
  • LOF は、PR-AUC 0.927 と飛び抜けて良く、上位 20 点のうち 16 点が異常で、すぐ外の異常も 9 点見つけました。距離そのものではなく、近所の点と比べてまばらかどうかを見るので、密なまとまりの近くでは小さなずれも浮いていると判定できるのです
  • Isolation Forest は、ROC-AUC 0.973 と全体の順位付けは良いものの、すぐ外の異常は 2 点しか見つけられませんでした。距離を計算しないので、データが大きくても速く動くのが長所です

ROC-AUC では 0.956〜0.998 と差が小さく見える方法どうしでも、PR-AUC では 0.485〜0.927 と大きな差があります。上位から調べていくときに、どれだけ効率よく本物の異常に当たるかは、PR-AUC のほうがよく表します。

万能な方法はない

この例では LOF が最も良い結果でしたが、いつもそうとは限りません。LOF は、全部の点どうしの距離を計算するので、データが大きいと重くなります。また、何次元もあるデータでは、どの点どうしの距離も同じくらいになってしまい、距離に基づく方法はうまく働かなくなります(次元の呪い)。データの形・大きさ・次元に合わせて、方法を選ぶ必要があります。

多くのセンサーの連動の崩れ

主成分分析で圧縮して戻す

工場の設備には、温度・圧力・振動・電流など、何十ものセンサーが付いています。これらの値は、「負荷」「気温」「回転数」のような少数の状態に合わせて、互いに連動して動きます。故障すると、この連動が崩れます。1 つのセンサーの値だけを見ると普通の範囲でも、ほかのセンサーとのつじつまが合わなくなるのです。

この崩れを見つけるには、クラスタリングと次元削減のしくみ で説明した 主成分分析(PCA)を使います。正常なときのデータで、20 個のセンサーの値を少数の主成分に圧縮して戻す方法を学びます。正常なデータは連動しているので、圧縮して戻してもほとんど元に戻りますが、連動が崩れたデータは、うまく戻りません。この戻しきれなかった分を 再構成誤差と呼び、異常の点数にします。

20 個のセンサーが 3 つの状態に合わせて動くデータで、1 つのセンサーだけが少しずれる故障を、センサーごとの z スコアと、再構成誤差で見つけられるかを比べます。しきい値は、学習に使っていない正常なデータの 1% だけが超える値にします(誤報の率を 1% にする)。

import numpy as np

# 20 個のセンサーの値。工場の 3 つの状態(負荷・気温・回転数)で、互いに連動して動く
rng = np.random.default_rng(0)
W = rng.standard_normal((3, 20))
def normal_data(n):
    return rng.standard_normal((n, 3)) @ W + 0.3 * rng.standard_normal((n, 20))
train, valid = normal_data(2000), normal_data(1000)        # 正常なときの記録だけで学習する
test_normal = normal_data(500)
test_fault = normal_data(500)
broken = rng.integers(20, size=500)                         # 故障:1 個のセンサーだけがずれる
test_fault[np.arange(500), broken] += rng.choice([-1, 1], 500) * 1.5

mu, sd = train.mean(axis=0), train.std(axis=0)
def score_z(X):                                             # センサーごとに見て、いちばん外れた値
    return np.abs((X - mu) / sd).max(axis=1)

# 主成分分析で 3 次元に圧縮して戻し、戻しきれなかった分(再構成誤差)を見る
_, _, Vt = np.linalg.svd(train - mu, full_matrices=False)
P = Vt[:3]
def score_pca(X):
    Z = (X - mu) @ P.T                                      # 圧縮
    return np.linalg.norm((X - mu) - Z @ P, axis=1)         # 戻したものとの差

print(round(float(np.abs(test_fault[np.arange(500), broken] - mu[broken]).mean() / sd[broken].mean()), 2))
for name, f in [("センサーごと", score_z), ("再構成誤差", score_pca)]:
    th = np.quantile(f(valid), 0.99)                        # 正常なデータの 1% だけが超える値をしきい値に
    print(name, round(float((f(test_normal) > th).mean()), 3), round(float((f(test_fault) > th).mean()), 3))

実行結果:

1.22
センサーごと 0.008 0.04
再構成誤差 0.01 0.674

1 行目は、故障したセンサーの値が平均から何標準偏差離れているかの平均で、1.22 です。故障したセンサーの値そのものは、正常な範囲に収まっています。

2・3 行目は、左から、正常なデータを誤って異常とした割合(誤報)と、故障を見つけた割合です。誤報はどちらも 1% 前後にそろえています。センサーを 1 つずつ見る方法では、故障を 4% しか見つけられませんでしたが、再構成誤差では 67.4% を見つけられました。値の大きさではなく、ほかのセンサーとのつじつまの崩れを見ているからです。

オートエンコーダー

主成分分析は「直線的な圧縮」しかできません。センサーの連動が曲がった関係のときは、ニューラルネットワークで圧縮して戻す オートエンコーダー(画像生成 AI のしくみ で説明したもの)を使います。正常な画像だけで学習したオートエンコーダーに、傷のある製品の画像を入れると、傷の部分がうまく戻らないので、その場所も分かります。外観検査では、正常な画像の特徴を学習済みの画像認識モデルで取り出し、そこからの k 近傍距離を測る方法(2022 年の PatchCore など)もよく使われています。

しきい値と運用

異常検知のモデルが出すのは「異常らしさ」の点数で、どこから異常とみなすかの しきい値は、人が決めます。しきい値を下げれば見逃しは減りますが、誤報が増えます。誤報が多すぎると、担当者が警報を無視するようになってしまうので、「1 日に何件までなら確認できるか」から逆算して決めることもよくあります。

また、正常なデータの姿は、季節・設備の交換・利用者の増加などで、少しずつ変わっていきます(ドリフト)。昔の正常なデータで学習したモデルは、新しい正常なデータを異常と判定し始めるので、定期的に学び直す必要があります。時間とともに変わるデータの扱い方は、時系列予測のしくみ も参考になります。予測した値と実際の値のずれを異常の点数にする方法も、よく使われます。

異常検知の手法の流れ

この記事で紹介した手法を年代順に並べます。

年手法・出来事ポイント
1936マハラノビス距離変数どうしの相関を考えた距離
1947ホテリングの T²多変量のデータの品質管理
1950グラブス検定正規分布を前提に、外れ値かどうかを検定する
1974ハンペルの頑健統計中央値・MAD など、外れ値に引っ張られにくい統計
2000LOF近所の点と比べた局所的な密度で外れ値を測る
20011 クラス SVM正常なデータを囲む境界を学ぶ
2008Isolation Forestランダムに切って、孤立させやすさで測る。距離を計算しないので速い
2014オートエンコーダーによる異常検知の広がり再構成誤差で、曲がった関係の崩れも見つける
2019MVTec AD製品の外観検査の異常検知のデータセット
2022PatchCore学習済みの画像の特徴からの近傍距離で、傷の場所を見つける

まとめ

  • 異常の例はほとんど集まらないので、正常なデータの姿を学び、そこから外れたものを異常とみなす
  • 1 つの値なら z スコアで外れ値を見つけられる。異常が多く混じるときは、中央値と MAD を使う頑健な方法が引っ張られにくい
  • 1 つずつ見ると普通でも、組み合わせがおかしいデータは、相関を考えたマハラノビス距離で見つけられる
  • 正常なデータがいくつものまとまりに分かれていると、全体の中心からの距離は役に立たない。k 近傍距離・LOF・Isolation Forest などを、データに合わせて選ぶ
  • 異常は少ないので、正解率ではなく PR-AUC などで評価する
  • 多くのセンサーが連動するデータでは、主成分分析やオートエンコーダーの再構成誤差で、連動の崩れを見つけられる
  • しきい値は誤報と見逃しの釣り合いで決め、正常なデータの変化に合わせて学び直す

参考

  • P. C. Mahalanobis, "On the generalized distance in statistics", Proceedings of the National Institute of Sciences of India 1936
  • H. Hotelling, "Multivariate Quality Control", Techniques of Statistical Analysis 1947
  • F. E. Grubbs, "Sample Criteria for Testing Outlying Observations", The Annals of Mathematical Statistics 1950
  • F. R. Hampel, "The Influence Curve and Its Role in Robust Estimation", Journal of the American Statistical Association 1974
  • M. M. Breunig et al., "LOF: Identifying Density-Based Local Outliers", SIGMOD 2000
  • B. Schölkopf et al., "Estimating the Support of a High-Dimensional Distribution", Neural Computation 2001
  • F. T. Liu, K. M. Ting and Z.-H. Zhou, "Isolation Forest", ICDM 2008
  • V. Chandola, A. Banerjee and V. Kumar, "Anomaly Detection: A Survey", ACM Computing Surveys 2009
  • M. Sakurada and T. Yairi, "Anomaly Detection Using Autoencoders with Nonlinear Dimensionality Reduction", MLSDA 2014
  • P. Bergmann et al., "MVTec AD — A Comprehensive Real-World Dataset for Unsupervised Anomaly Detection", CVPR 2019
  • K. Roth et al., "Towards Total Recall in Industrial Anomaly Detection", CVPR 2022(arXiv:2106.08265)