製造業における不良品検出、インフラの故障予知、サイバーセキュリティにおける侵入検知など、異常検知は多くの実務で必要とされるタスクです。しかし、現実世界では「異常データが極端に少ない」という問題があり、教師あり学習の適用が困難なケースがほとんどです。
オートエンコーダ（Autoencoder, AE）による異常検知は、「正常データだけで学習し、正常パターンを再構成できるモデルを構築する」というアプローチをとります。異常データは学習時に見ていないため再構成がうまくいかず、再構成誤差が大きくなることを利用して異常を検出します。
本記事の内容
- AEベース異常検知の原理と再構成誤差
- 閾値設定の手法（統計的アプローチ・極値理論）
- ボトルネック次元の選択指針
- 畳み込みAEによる画像異常検知の概要
- LSTMオートエンコーダによる時系列異常検知の概要
- アンサンブルAEの概要
- Pythonでの合成製造データによるAE異常検知の実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
AEベース異常検知の原理
オートエンコーダの復習
オートエンコーダは、エンコーダ $E_\phi$ とデコーダ $D_\theta$ から構成されるニューラルネットワークです。
$$ \begin{align} \bm{z} &= E_\phi(\bm{x}) \quad \text{(エンコーダ: 入力を潜在表現に変換)} \\ \hat{\bm{x}} &= D_\theta(\bm{z}) \quad \text{(デコーダ: 潜在表現から入力を再構成)} \end{align} $$
ここで $\bm{x} \in \mathbb{R}^d$ は入力データ、$\bm{z} \in \mathbb{R}^k$ は潜在表現（$k \ll d$）、$\hat{\bm{x}} \in \mathbb{R}^d$ は再構成データです。
学習の目的関数は再構成誤差の最小化です。
$$ \min_{\phi, \theta} \frac{1}{N}\sum_{i=1}^{N} \|\bm{x}_i – D_\theta(E_\phi(\bm{x}_i))\|^2 $$
異常検知への応用
AEベース異常検知のアイデアは非常にシンプルです。
- 学習フェーズ: 正常データのみでオートエンコーダを学習する
- 推論フェーズ: 新しいデータ $\bm{x}$ を入力し、再構成 $\hat{\bm{x}}$ を得る
- 異常判定: 再構成誤差 $e(\bm{x})$ が閾値を超えたら異常と判定する
この手法が機能する理由は、ボトルネック構造にあります。潜在空間の次元 $k$ が入力次元 $d$ より小さいため、AEは入力の全ての情報を保持できません。学習データ（正常データ）に多いパターンの再構成が優先的に学習されるため、未知の異常パターンは再構成精度が低くなります。
数学的な理解
正常データの分布を $p_{\text{normal}}(\bm{x})$ とすると、AEの学習は近似的に以下を実現します。
$$ D_\theta(E_\phi(\bm{x})) \approx \mathbb{E}[\bm{x} | \bm{z}], \quad \bm{z} = E_\phi(\bm{x}) $$
すなわち、潜在表現 $\bm{z}$ が与えられたときの入力の条件付き期待値を出力します。正常データは潜在空間上で適切にエンコードされるため、この条件付き期待値が入力自身に近くなります。一方、異常データは正常データの潜在空間の「外」にマッピングされるか、正常データとは異なる領域にマッピングされるため、再構成精度が低下します。
再構成誤差による異常スコア
基本的な再構成誤差
最も基本的な異常スコアは、入力と再構成の間のL2ノルム（二乗誤差）です。
$$ e(\bm{x}) = \|\bm{x} – D_\theta(E_\phi(\bm{x}))\|_2^2 = \sum_{j=1}^{d} (x_j – \hat{x}_j)^2 $$
正規化された再構成誤差
各特徴量のスケールが異なる場合、正規化が必要です。学習データでの各特徴量の再構成誤差の平均 $\mu_j$ と標準偏差 $\sigma_j$ を用いて、
$$ e_{\text{norm}}(\bm{x}) = \sum_{j=1}^{d} \left(\frac{(x_j – \hat{x}_j)^2 – \mu_j}{\sigma_j}\right)^2 $$
ここで
$$ \mu_j = \frac{1}{N}\sum_{i=1}^{N}(x_{ij} – \hat{x}_{ij})^2, \quad \sigma_j^2 = \frac{1}{N}\sum_{i=1}^{N}\left((x_{ij} – \hat{x}_{ij})^2 – \mu_j\right)^2 $$
特徴量ごとの再構成誤差
異常の原因を特定するために、特徴量ごとの再構成誤差も有用です。
$$ e_j(\bm{x}) = (x_j – \hat{x}_j)^2 $$
$e_j(\bm{x})$ が大きい特徴量 $j$ が、異常の原因である可能性が高いと解釈できます。
閾値設定の手法
統計的アプローチ
正常データでの再構成誤差の分布を用いて閾値を設定する最も単純な方法です。
パーセンタイルベース: 学習データの再構成誤差の上位 $\alpha$% パーセンタイルを閾値とします。
$$ \tau = Q_{1-\alpha}(\{e(\bm{x}_i)\}_{i=1}^{N}) $$
例えば $\alpha = 0.01$ なら、正常データの99パーセンタイルを閾値とし、正常データの1%を擬似的な異常として許容します。
ガウス仮定: 再構成誤差が正規分布に従うと仮定する場合、
$$ e(\bm{x}) \sim \mathcal{N}(\mu_e, \sigma_e^2) $$
閾値は $\tau = \mu_e + k\sigma_e$ とします。$k = 3$ なら約0.3%の正常データが異常と判定されます。
極値理論（EVT）による閾値設定
再構成誤差が正規分布に従わない場合、極値理論（Extreme Value Theory）に基づくアプローチが有効です。
一般化パレート分布（GPD）
EVTのPOT（Peaks Over Threshold）法では、十分に高い閾値 $u$ を超えた部分の分布が一般化パレート分布（GPD）に従うことを利用します。
$$ P(e > u + y | e > u) \approx G_{\gamma, \sigma}(y) = 1 – \left(1 + \gamma\frac{y}{\sigma}\right)^{-1/\gamma} $$
ここで $\gamma$ は形状パラメータ、$\sigma$ はスケールパラメータです。
POT法の手順
- 初期閾値 $u$ を設定する（例: 再構成誤差の90パーセンタイル）
- $u$ を超える再構成誤差の超過量 $y_i = e(\bm{x}_i) – u$ を集める
- 超過量にGPDをフィットし、パラメータ $(\hat{\gamma}, \hat{\sigma})$ を推定する
- 所望の異常率 $q$ に対応する閾値 $\tau$ を計算する
最終的な閾値は次のように求められます。
$$ \tau = u + \frac{\hat{\sigma}}{\hat{\gamma}}\left[\left(\frac{n}{N_u} \cdot \frac{1}{q}\right)^{\hat{\gamma}} – 1\right] $$
ここで $n$ は全データ数、$N_u$ は閾値 $u$ を超えたデータ数、$q$ は目標とする偽陽性率です。
GPDパラメータの最尤推定
超過量 $y_1, y_2, \dots, y_{N_u}$ に対するGPDの対数尤度は、
$$ \ell(\gamma, \sigma) = -N_u \ln \sigma – \left(1 + \frac{1}{\gamma}\right)\sum_{i=1}^{N_u}\ln\left(1 + \gamma\frac{y_i}{\sigma}\right) $$
この対数尤度を最大化するパラメータ $(\hat{\gamma}, \hat{\sigma})$ を数値最適化で求めます。
偏微分を計算すると、
$$ \begin{align} \frac{\partial \ell}{\partial \sigma} &= -\frac{N_u}{\sigma} + \left(1 + \frac{1}{\gamma}\right)\sum_{i=1}^{N_u}\frac{\gamma y_i / \sigma^2}{1 + \gamma y_i / \sigma} = 0 \\ \frac{\partial \ell}{\partial \gamma} &= \frac{1}{\gamma^2}\sum_{i=1}^{N_u}\ln\left(1 + \gamma\frac{y_i}{\sigma}\right) – \left(1 + \frac{1}{\gamma}\right)\sum_{i=1}^{N_u}\frac{y_i / \sigma}{1 + \gamma y_i / \sigma} = 0 \end{align} $$
この連立方程式は解析的に解けないため、Newton-Raphson法やscipy.optimize.minimizeなどの数値最適化を用います。
ボトルネック次元の選択
次元選択の基本原理
ボトルネック次元 $k$ はAEベース異常検知の性能を大きく左右します。
- $k$ が大きすぎる: 恒等写像に近くなり、異常データも正確に再構成してしまう
- $k$ が小さすぎる: 正常データの再構成精度も低くなり、正常と異常の区別が困難になる
累積寄与率との関係
PCA（主成分分析）との類似性から、ボトルネック次元の選択指針を得ることができます。線形AEはPCAと等価であることが知られています。
データの共分散行列の固有値を $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_d$ とすると、累積寄与率は
$$ R(k) = \frac{\sum_{i=1}^{k} \lambda_i}{\sum_{i=1}^{d} \lambda_i} $$
$R(k)$ が90-95%程度になる $k$ を目安にボトルネック次元を選択します。ただし、非線形AEはPCAより表現力が高いため、実際にはこれより小さい $k$ でも十分な再構成精度を達成できることが多いです。
実用的なアプローチ
- 複数の $k$ でAEを学習し、検証データでの異常検知性能（AUC-ROCなど）を比較する
- 再構成誤差のエルボー法: $k$ を増やしていき、再構成誤差の減少が緩やかになる点を選ぶ
畳み込みAEによる画像異常検知
画像データの異常検知には、畳み込みオートエンコーダ（Convolutional Autoencoder, CAE）が有効です。
アーキテクチャ
エンコーダは畳み込み層とプーリング層で空間的な特徴を抽出し、デコーダは転置畳み込み（アップサンプリング）で画像を再構成します。
$$ \begin{align} \bm{z} &= \text{Conv}_{\text{enc}}(\bm{x}) \quad \text{(畳み込み + プーリング)} \\ \hat{\bm{x}} &= \text{ConvT}_{\text{dec}}(\bm{z}) \quad \text{(転置畳み込み + アップサンプリング)} \end{align} $$
画像異常検知の特徴
画像異常検知では、ピクセルレベルの再構成誤差マップを計算することで、異常の位置を特定できます。
$$ \bm{M}(\bm{x}) = |\bm{x} – \hat{\bm{x}}|^2 $$
この誤差マップ $\bm{M}(\bm{x})$ の各ピクセル値が高い領域が、異常が存在する位置を示します。
SSIM（Structural Similarity Index）ベースの異常スコア
ピクセル単位のL2誤差の代わりに、SSIM（Structural Similarity Index）を用いた異常スコアも有効です。
$$ \text{SSIM}(\bm{x}, \hat{\bm{x}}) = \frac{(2\mu_x\mu_{\hat{x}} + c_1)(2\sigma_{x\hat{x}} + c_2)}{(\mu_x^2 + \mu_{\hat{x}}^2 + c_1)(\sigma_x^2 + \sigma_{\hat{x}}^2 + c_2)} $$
ここで $\mu_x, \mu_{\hat{x}}$ は局所平均、$\sigma_x^2, \sigma_{\hat{x}}^2$ は局所分散、$\sigma_{x\hat{x}}$ は局所共分散、$c_1, c_2$ は安定化定数です。SSIMは人間の知覚に近い構造的な類似度を測ることができます。
異常スコアは $1 – \text{SSIM}$ として定義します。
LSTMオートエンコーダによる時系列異常検知
時系列データの異常検知には、LSTM（Long Short-Term Memory）オートエンコーダが有効です。
アーキテクチャの概要
LSTMオートエンコーダでは、エンコーダがLSTM層で時系列の情報を潜在ベクトル（最終隠れ状態）に圧縮し、デコーダがLSTM層で時系列を再構成します。
入力時系列 $\bm{x}_{1:T} = (\bm{x}_1, \bm{x}_2, \dots, \bm{x}_T)$ に対して、
$$ \begin{align} \bm{h}_T &= \text{LSTM}_{\text{enc}}(\bm{x}_{1:T}) \quad \text{(潜在表現)} \\ \hat{\bm{x}}_{1:T} &= \text{LSTM}_{\text{dec}}(\bm{h}_T) \quad \text{(再構成時系列)} \end{align} $$
時系列の異常スコア
時系列全体の再構成誤差は、
$$ e(\bm{x}_{1:T}) = \frac{1}{T}\sum_{t=1}^{T}\|\bm{x}_t – \hat{\bm{x}}_t\|^2 $$
各時刻の再構成誤差 $e_t = \|\bm{x}_t – \hat{\bm{x}}_t\|^2$ をモニタリングすることで、異常が発生した時刻を特定することもできます。
アンサンブルAE
基本的なアイデア
複数のAEを学習し、その異常スコアを集約することで、単一のAEよりもロバストな異常検知が可能になります。
$M$ 個のAE $\{(E_{\phi_m}, D_{\theta_m})\}_{m=1}^{M}$ を学習し、アンサンブル異常スコアを
$$ e_{\text{ensemble}}(\bm{x}) = \frac{1}{M}\sum_{m=1}^{M} e_m(\bm{x}) = \frac{1}{M}\sum_{m=1}^{M}\|\bm{x} – D_{\theta_m}(E_{\phi_m}(\bm{x}))\|^2 $$
として定義します。
多様性の確保
アンサンブルの効果を高めるためには、各AEの多様性を確保することが重要です。
- 異なるアーキテクチャ: ボトルネック次元、層数、活性化関数を変える
- 異なるサブセット: 学習データのサブセットを変える（Bagging的アプローチ）
- 異なる特徴量: 入力特徴量のサブセットを変える（Random Subspace法）
- 異なる初期化: 同じアーキテクチャでも重みの初期化を変える
Pythonでの実装
合成製造データによるAE異常検知
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import roc_auc_score, roc_curve, precision_recall_curve
# ===== 1. 合成製造データの生成 =====
np.random.seed(42)
def generate_manufacturing_data(n_normal=1000, n_anomaly=50, n_features=10):
"""
製造データを模した合成データを生成
正常データ: 相関のある多変量正規分布
異常データ: 一部の特徴量にシフトや異常値を含む
"""
# 正常データの共分散行列を生成（相関のある特徴量）
A = np.random.randn(n_features, n_features) * 0.3
cov = A @ A.T + np.eye(n_features) * 0.5
# 正常データの平均
mean_normal = np.zeros(n_features)
# 正常データの生成
X_normal = np.random.multivariate_normal(mean_normal, cov, n_normal)
# 異常データの生成（複数の異常パターン）
X_anomaly = np.random.multivariate_normal(mean_normal, cov, n_anomaly)
# 異常パターン1: 一部特徴量のシフト（40%）
n_shift = int(n_anomaly * 0.4)
for i in range(n_shift):
shift_features = np.random.choice(n_features, 3, replace=False)
X_anomaly[i, shift_features] += np.random.uniform(3, 5, 3) * np.random.choice([-1, 1], 3)
# 異常パターン2: 外れ値（30%）
n_outlier = int(n_anomaly * 0.3)
for i in range(n_shift, n_shift + n_outlier):
outlier_features = np.random.choice(n_features, 2, replace=False)
X_anomaly[i, outlier_features] = np.random.uniform(5, 8, 2) * np.random.choice([-1, 1], 2)
# 異常パターン3: 相関構造の破壊（30%）
n_corr_break = n_anomaly - n_shift - n_outlier
for i in range(n_shift + n_outlier, n_anomaly):
X_anomaly[i] = np.random.randn(n_features) * 3
# ラベル
y = np.array([0] * n_normal + [1] * n_anomaly)
return np.vstack([X_normal, X_anomaly]), y
X, y_true = generate_manufacturing_data(n_normal=1000, n_anomaly=50, n_features=10)
print(f"データ形状: {X.shape}")
print(f"正常: {np.sum(y_true == 0)}, 異常: {np.sum(y_true == 1)}")
# 標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
# 学習/テスト分割（正常データのみで学習）
X_train_normal = X_scaled[y_true == 0][:800] # 正常データの80%で学習
X_val_normal = X_scaled[y_true == 0][800:] # 正常データの20%で検証
X_test = X_scaled # 全データでテスト
y_test = y_true
# ===== 2. オートエンコーダの実装（NumPyベース） =====
class SimpleAutoencoder:
"""NumPyベースの全結合オートエンコーダ"""
def __init__(self, input_dim, hidden_dims=[64, 32, 16, 32, 64],
learning_rate=0.001):
"""
input_dim: 入力次元
hidden_dims: 隠れ層の次元リスト（エンコーダ→ボトルネック→デコーダ）
"""
self.lr = learning_rate
self.layers = []
dims = [input_dim] + hidden_dims + [input_dim]
# 重みの初期化（Xavier初期化）
self.weights = []
self.biases = []
for i in range(len(dims) - 1):
scale = np.sqrt(2.0 / (dims[i] + dims[i+1]))
W = np.random.randn(dims[i], dims[i+1]) * scale
b = np.zeros(dims[i+1])
self.weights.append(W)
self.biases.append(b)
self.n_layers = len(self.weights)
def relu(self, x):
return np.maximum(0, x)
def relu_deriv(self, x):
return (x > 0).astype(float)
def forward(self, X):
"""順伝播"""
self.activations = [X]
self.pre_activations = []
h = X
for i in range(self.n_layers):
z = h @ self.weights[i] + self.biases[i]
self.pre_activations.append(z)
if i < self.n_layers - 1:
h = self.relu(z) # 中間層はReLU
else:
h = z # 出力層は線形
self.activations.append(h)
return h
def backward(self, X, X_hat):
"""逆伝播"""
m = X.shape[0]
# 出力層の勾配（MSE損失の勾配）
delta = 2 * (X_hat - X) / m
grad_weights = []
grad_biases = []
for i in range(self.n_layers - 1, -1, -1):
grad_W = self.activations[i].T @ delta
grad_b = np.sum(delta, axis=0)
grad_weights.insert(0, grad_W)
grad_biases.insert(0, grad_b)
if i > 0:
delta = (delta @ self.weights[i].T) * self.relu_deriv(
self.pre_activations[i-1]
)
return grad_weights, grad_biases
def update(self, grad_weights, grad_biases):
"""パラメータ更新（SGD）"""
for i in range(self.n_layers):
self.weights[i] -= self.lr * grad_weights[i]
self.biases[i] -= self.lr * grad_biases[i]
def fit(self, X_train, epochs=200, batch_size=64, verbose=True):
"""学習"""
n = X_train.shape[0]
losses = []
for epoch in range(epochs):
# ミニバッチでシャッフル
indices = np.random.permutation(n)
epoch_loss = 0
for start in range(0, n, batch_size):
end = min(start + batch_size, n)
X_batch = X_train[indices[start:end]]
# 順伝播
X_hat = self.forward(X_batch)
# 損失計算
loss = np.mean((X_batch - X_hat) ** 2)
epoch_loss += loss * (end - start)
# 逆伝播
grad_W, grad_b = self.backward(X_batch, X_hat)
# 更新
self.update(grad_W, grad_b)
epoch_loss /= n
losses.append(epoch_loss)
if verbose and (epoch + 1) % 50 == 0:
print(f"Epoch {epoch+1}/{epochs}, Loss: {epoch_loss:.6f}")
return losses
def reconstruct(self, X):
"""再構成"""
return self.forward(X)
def reconstruction_error(self, X):
"""再構成誤差を計算"""
X_hat = self.reconstruct(X)
return np.mean((X - X_hat) ** 2, axis=1)
print("SimpleAutoencoder クラスの定義完了")
# ===== 3. モデルの学習 =====
input_dim = X_scaled.shape[1]
ae = SimpleAutoencoder(
input_dim=input_dim,
hidden_dims=[64, 32, 8, 32, 64], # ボトルネック次元=8
learning_rate=0.005
)
losses = ae.fit(X_train_normal, epochs=300, batch_size=64, verbose=True)
# 学習曲線の可視化
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(losses, linewidth=2)
ax.set_xlabel('Epoch', fontsize=12)
ax.set_ylabel('MSE Loss', fontsize=12)
ax.set_title('Autoencoder Training Loss', fontsize=14)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# ===== 4. 異常検知と閾値設定 =====
# テストデータの再構成誤差
errors_test = ae.reconstruction_error(X_test)
errors_train = ae.reconstruction_error(X_train_normal)
errors_val = ae.reconstruction_error(X_val_normal)
# 方法1: パーセンタイルベースの閾値
threshold_percentile = np.percentile(errors_train, 99)
print(f"閾値（99パーセンタイル）: {threshold_percentile:.6f}")
# 方法2: 平均 + k*標準偏差
mu_e = np.mean(errors_train)
sigma_e = np.std(errors_train)
threshold_gaussian = mu_e + 3 * sigma_e
print(f"閾値（平均+3σ）: {threshold_gaussian:.6f}")
# 方法3: EVT（POT法）による閾値
def fit_gpd_pot(errors, initial_threshold_percentile=90, target_fpr=0.01):
"""
POT法による閾値設定
"""
from scipy.optimize import minimize
from scipy.stats import genpareto
# 初期閾値
u = np.percentile(errors, initial_threshold_percentile)
exceedances = errors[errors > u] - u
if len(exceedances) < 10:
print("警告: 超過データが少なすぎます。パーセンタイル法にフォールバック")
return np.percentile(errors, (1 - target_fpr) * 100)
# GPDのフィット
shape, loc, scale = genpareto.fit(exceedances, floc=0)
# 閾値の計算
n = len(errors)
n_u = len(exceedances)
tau = u + (scale / shape) * ((n / n_u * (1 / target_fpr)) ** shape - 1)
return tau, shape, scale, u
threshold_evt, gamma_hat, sigma_hat, u_pot = fit_gpd_pot(errors_train)
print(f"閾値（EVT-POT）: {threshold_evt:.6f}")
print(f" GPDパラメータ: γ={gamma_hat:.4f}, σ={sigma_hat:.4f}, u={u_pot:.6f}")
# 各閾値での判定結果
thresholds = {
'99パーセンタイル': threshold_percentile,
'平均+3σ': threshold_gaussian,
'EVT-POT': threshold_evt
}
print("\n各閾値での異常検知性能:")
print(f"{'方法':<18} {'閾値':>10} {'検出数':>8} {'真の異常数':>10} {'精度':>8}")
print("-" * 58)
for name, thresh in thresholds.items():
pred = (errors_test > thresh).astype(int)
n_detected = np.sum(pred)
n_true_detected = np.sum(pred[y_test == 1])
n_anomaly = np.sum(y_test == 1)
precision = n_true_detected / max(n_detected, 1)
print(f"{name:<18} {thresh:>10.4f} {n_detected:>8d} {n_true_detected:>4d}/{n_anomaly:<5d} {precision:>8.4f}")
# ===== 5. 可視化と評価 =====
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (a) 再構成誤差の分布
axes[0, 0].hist(errors_test[y_test == 0], bins=50, alpha=0.7, label='Normal',
density=True, color='blue')
axes[0, 0].hist(errors_test[y_test == 1], bins=30, alpha=0.7, label='Anomaly',
density=True, color='red')
for name, thresh in thresholds.items():
axes[0, 0].axvline(x=thresh, linestyle='--', label=f'{name}: {thresh:.3f}')
axes[0, 0].set_xlabel('Reconstruction Error')
axes[0, 0].set_ylabel('Density')
axes[0, 0].set_title('Reconstruction Error Distribution')
axes[0, 0].legend(fontsize=8)
axes[0, 0].grid(True, alpha=0.3)
# (b) ROC曲線
fpr, tpr, _ = roc_curve(y_test, errors_test)
auc = roc_auc_score(y_test, errors_test)
axes[0, 1].plot(fpr, tpr, 'b-', lw=2, label=f'AE (AUC={auc:.4f})')
axes[0, 1].plot([0, 1], [0, 1], 'k--', lw=1, alpha=0.5)
axes[0, 1].set_xlabel('False Positive Rate')
axes[0, 1].set_ylabel('True Positive Rate')
axes[0, 1].set_title('ROC Curve')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# (c) Precision-Recall曲線
precision, recall, _ = precision_recall_curve(y_test, errors_test)
axes[1, 0].plot(recall, precision, 'g-', lw=2)
axes[1, 0].set_xlabel('Recall')
axes[1, 0].set_ylabel('Precision')
axes[1, 0].set_title('Precision-Recall Curve')
axes[1, 0].grid(True, alpha=0.3)
# (d) 特徴量ごとの再構成誤差（上位異常データ）
X_hat_test = ae.reconstruct(X_test)
feature_errors = (X_test - X_hat_test) ** 2
# 異常スコアが高い上位5点
top_anomaly_idx = np.argsort(errors_test)[-5:]
for idx in top_anomaly_idx:
label = 'Anomaly' if y_test[idx] == 1 else 'Normal'
axes[1, 1].plot(range(input_dim), feature_errors[idx],
marker='o', label=f'#{idx} ({label})', alpha=0.7)
axes[1, 1].set_xlabel('Feature Index')
axes[1, 1].set_ylabel('Reconstruction Error')
axes[1, 1].set_title('Per-Feature Error (Top 5 Anomalies)')
axes[1, 1].legend(fontsize=8)
axes[1, 1].grid(True, alpha=0.3)
plt.suptitle('Autoencoder Anomaly Detection Results', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
print(f"\nAUC-ROC: {auc:.4f}")
# ===== 6. ボトルネック次元の影響分析 =====
bottleneck_dims = [2, 4, 8, 16, 32]
auc_by_dim = []
for bdim in bottleneck_dims:
ae_test = SimpleAutoencoder(
input_dim=input_dim,
hidden_dims=[64, 32, bdim, 32, 64],
learning_rate=0.005
)
ae_test.fit(X_train_normal, epochs=200, batch_size=64, verbose=False)
errors = ae_test.reconstruction_error(X_test)
auc_val = roc_auc_score(y_test, errors)
train_error = np.mean(ae_test.reconstruction_error(X_train_normal))
auc_by_dim.append({
'dim': bdim,
'auc': auc_val,
'train_error': train_error
})
print(f"ボトルネック次元={bdim}: AUC={auc_val:.4f}, 学習誤差={train_error:.6f}")
# 可視化
fig, ax1 = plt.subplots(figsize=(8, 5))
dims = [r['dim'] for r in auc_by_dim]
aucs = [r['auc'] for r in auc_by_dim]
train_errors = [r['train_error'] for r in auc_by_dim]
ax1.plot(dims, aucs, 'bo-', lw=2, label='AUC-ROC')
ax1.set_xlabel('Bottleneck Dimension', fontsize=12)
ax1.set_ylabel('AUC-ROC', fontsize=12, color='blue')
ax1.tick_params(axis='y', labelcolor='blue')
ax2 = ax1.twinx()
ax2.plot(dims, train_errors, 'rs--', lw=2, label='Train Error')
ax2.set_ylabel('Train Reconstruction Error', fontsize=12, color='red')
ax2.tick_params(axis='y', labelcolor='red')
ax1.set_title('Effect of Bottleneck Dimension', fontsize=14)
ax1.grid(True, alpha=0.3)
fig.legend(loc='upper center', bbox_to_anchor=(0.5, 0.92), ncol=2)
plt.tight_layout()
plt.show()
# ===== 7. アンサンブルAEの実装と評価 =====
n_ensemble = 5
ensemble_errors = np.zeros((X_test.shape[0], n_ensemble))
for m in range(n_ensemble):
# 異なるアーキテクチャ（ボトルネック次元をランダムに変える）
bdim = np.random.choice([4, 8, 12, 16])
ae_m = SimpleAutoencoder(
input_dim=input_dim,
hidden_dims=[64, 32, bdim, 32, 64],
learning_rate=0.005
)
# ブートストラップサンプリング
bootstrap_idx = np.random.choice(len(X_train_normal),
size=len(X_train_normal), replace=True)
X_train_bootstrap = X_train_normal[bootstrap_idx]
ae_m.fit(X_train_bootstrap, epochs=200, batch_size=64, verbose=False)
ensemble_errors[:, m] = ae_m.reconstruction_error(X_test)
print(f"アンサンブルモデル {m+1}/{n_ensemble} (bottleneck={bdim}): 完了")
# アンサンブルスコア（平均）
ensemble_score = np.mean(ensemble_errors, axis=1)
# 単一AEとアンサンブルの比較
auc_single = roc_auc_score(y_test, errors_test)
auc_ensemble = roc_auc_score(y_test, ensemble_score)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# ROC曲線の比較
fpr_s, tpr_s, _ = roc_curve(y_test, errors_test)
fpr_e, tpr_e, _ = roc_curve(y_test, ensemble_score)
axes[0].plot(fpr_s, tpr_s, 'b-', lw=2, label=f'Single AE (AUC={auc_single:.4f})')
axes[0].plot(fpr_e, tpr_e, 'r-', lw=2, label=f'Ensemble AE (AUC={auc_ensemble:.4f})')
axes[0].plot([0, 1], [0, 1], 'k--', lw=1, alpha=0.5)
axes[0].set_xlabel('False Positive Rate')
axes[0].set_ylabel('True Positive Rate')
axes[0].set_title('Single AE vs Ensemble AE')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# スコア分布の比較
axes[1].hist(ensemble_score[y_test == 0], bins=50, alpha=0.7, label='Normal',
density=True, color='blue')
axes[1].hist(ensemble_score[y_test == 1], bins=30, alpha=0.7, label='Anomaly',
density=True, color='red')
axes[1].set_xlabel('Ensemble Anomaly Score')
axes[1].set_ylabel('Density')
axes[1].set_title('Ensemble Score Distribution')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"単一AE AUC: {auc_single:.4f}")
print(f"アンサンブルAE AUC: {auc_ensemble:.4f}")
まとめ
本記事では、オートエンコーダによる異常検知の理論と実装について解説しました。
- 基本原理: 正常データのみで学習したAEは、異常データの再構成誤差が大きくなることを利用して異常を検出する
- 再構成誤差: $e(\bm{x}) = \|\bm{x} – D_\theta(E_\phi(\bm{x}))\|^2$ がシンプルかつ効果的な異常スコアになる
- 閾値設定: パーセンタイル法、ガウス仮定法、極値理論（EVT）のPOT法の3つのアプローチを解説した
- ボトルネック次元: 小さすぎても大きすぎても性能が低下するため、適切な選択が重要
- 発展的手法: 畳み込みAE（画像）、LSTMオートエンコーダ（時系列）、アンサンブルAE（ロバスト性向上）
- Python実装: 合成製造データでAE異常検知を実装し、複数の閾値設定法と精度評価を実施した
次のステップとして、以下の記事も参考にしてください。