異常偵測：Isolation Forest、重建誤差與 GANomaly

本文使用的異常分數 = -score_samples(x)
decision_function(x) = score_samples(x) - offset_

contamination="auto" 時 offset_ = -0.5
本文不使用 predict 的預設分界，另選正常校準分數的閾值。

z = Encoder(x)
x_reconstructed = Decoder(z)
score_AE(x) = mean((x_reconstructed - x)^2)

這裡的 x 是按正常訓練資料 mean/std 縮放後的輸入。

z = Encoder_1(x)
x_reconstructed = Decoder(z)
z_reconstructed = Encoder_2(x_reconstructed)

本文概念模型的 generator objective：
feature_MSE + 50 * contextual_mean_L1 + latent_MSE

本文測試分數：mean((z - z_reconstructed)^2)

threshold = np.quantile(normal_calibration_scores, 0.95, method="higher")
flagged = test_scores > threshold

numpy==1.26.3
scikit-learn==1.5.1
torch==2.3.1

python3 -m venv .venv
source .venv/bin/activate
python -m pip install -r requirements.txt
python anomaly_scores_and_thresholds.py > stdout.json

"""Actual normal-only fits and separate normal calibration on a hand-set 2D toy.

EDE is a small encoder-decoder-encoder / feature-matching teaching model, not
a replication of GANomaly's image architecture or benchmark. No external data.
"""
import json
import platform

import numpy as np
import sklearn
import torch
from sklearn.ensemble import IsolationForest
from torch import nn


UPDATES = 500
SEED = 7


def normal_points(count, phase):
    t = -1 + 2 * (np.arange(count) + phase) / count
    disturbance = 0.02 * (np.arange(count) % 5 - 2)
    return np.column_stack((t, 0.5 * t + 0.15 * np.sin(3 * t) + disturbance))


class Autoencoder(nn.Module):
    def __init__(self):
        super().__init__()
        self.encoder = nn.Sequential(nn.Linear(2, 8), nn.Tanh(), nn.Linear(8, 1))
        self.decoder = nn.Sequential(nn.Linear(1, 8), nn.Tanh(), nn.Linear(8, 2))

    def forward(self, x):
        return self.decoder(self.encoder(x))


class EDE(nn.Module):
    def __init__(self):
        super().__init__()
        self.encoder1 = nn.Sequential(nn.Linear(2, 8), nn.Tanh(), nn.Linear(8, 1))
        self.decoder = nn.Sequential(nn.Linear(1, 8), nn.Tanh(), nn.Linear(8, 2))
        self.encoder2 = nn.Sequential(nn.Linear(2, 8), nn.Tanh(), nn.Linear(8, 1))

    def forward(self, x):
        z = self.encoder1(x)
        reconstructed = self.decoder(z)
        return reconstructed, z, self.encoder2(reconstructed)


class Discriminator(nn.Module):
    def __init__(self):
        super().__init__()
        self.features = nn.Sequential(nn.Linear(2, 8), nn.Tanh(), nn.Linear(8, 4), nn.Tanh())
        self.head = nn.Linear(4, 1)

    def forward(self, x):
        return self.head(self.features(x)).squeeze(-1)


def ede_losses(generator, discriminator, x):
    reconstructed, z, z_reconstructed = generator(x)
    feature = nn.functional.mse_loss(discriminator.features(reconstructed),
                                     discriminator.features(x).detach())
    contextual = nn.functional.l1_loss(reconstructed, x)
    latent = nn.functional.mse_loss(z_reconstructed, z)
    return feature, contextual, latent


def loss_record(generator, discriminator, x):
    with torch.no_grad():
        a, c, e = ede_losses(generator, discriminator, x)
    return {"feature_MSE": float(a), "contextual_L1_mean": float(c),
            "latent_MSE": float(e), "weighted_generator_objective": float(a + 50*c + e)}


def fit_networks(x):
    torch.manual_seed(SEED)
    autoencoder = Autoencoder().double()
    ae_optimizer = torch.optim.Adam(autoencoder.parameters(), lr=0.01)
    ae_initial = float(nn.functional.mse_loss(autoencoder(x), x).detach())
    for _ in range(UPDATES):
        ae_optimizer.zero_grad()
        loss = nn.functional.mse_loss(autoencoder(x), x)
        loss.backward()
        ae_optimizer.step()
    ae_final = float(nn.functional.mse_loss(autoencoder(x), x).detach())

    torch.manual_seed(SEED)
    generator, discriminator = EDE().double(), Discriminator().double()
    g_optimizer = torch.optim.Adam(generator.parameters(), lr=0.002, betas=(0.5, 0.999))
    d_optimizer = torch.optim.Adam(discriminator.parameters(), lr=0.002, betas=(0.5, 0.999))
    initial = loss_record(generator, discriminator, x)
    for _ in range(UPDATES):
        d_optimizer.zero_grad()
        with torch.no_grad():
            reconstructed = generator(x)[0]
        real_loss = nn.functional.binary_cross_entropy_with_logits(
            discriminator(x), torch.ones(len(x), dtype=torch.float64))
        fake_loss = nn.functional.binary_cross_entropy_with_logits(
            discriminator(reconstructed), torch.zeros(len(x), dtype=torch.float64))
        (0.5 * (real_loss + fake_loss)).backward()
        d_optimizer.step()
        # D weights stay fixed in this phase; gradients still pass through D to G.
        for parameter in discriminator.parameters():
            parameter.requires_grad_(False)
        g_optimizer.zero_grad()
        feature, contextual, latent = ede_losses(generator, discriminator, x)
        (feature + 50 * contextual + latent).backward()
        g_optimizer.step()
        for parameter in discriminator.parameters():
            parameter.requires_grad_(True)
    final = loss_record(generator, discriminator, x)
    autoencoder.eval()
    generator.eval()
    discriminator.eval()
    return autoencoder, generator, {"AE_initial_training_MSE": ae_initial,
        "AE_final_training_MSE": ae_final, "EDE_initial_training_losses": initial,
        "EDE_final_training_losses": final, "AE_updates": UPDATES,
        "generator_updates": UPDATES, "discriminator_updates": UPDATES,
        "parameters": {"autoencoder": sum(p.numel() for p in autoencoder.parameters()),
                       "EDE_generator": sum(p.numel() for p in generator.parameters()),
                       "discriminator": sum(p.numel() for p in discriminator.parameters())}}


def evaluate_scores(name, calibration, scores, labels):
    # Freeze the model-specific cutoff using normal calibration points only.
    threshold = float(np.quantile(calibration, 0.95, method="higher"))
    flagged = scores > threshold
    positive = labels == 1
    tp = int(np.sum(flagged & positive))
    fp = int(np.sum(flagged & ~positive))
    fn = int(np.sum(~flagged & positive))
    tn = int(np.sum(~flagged & ~positive))
    return {"model": name, "normal_calibration_scores": calibration.tolist(),
        "threshold": threshold, "threshold_rule": "np.quantile(normal_calibration, .95, method='higher'); flag score > threshold",
        "normal_calibration_flagged": int(np.sum(calibration > threshold)),
        "test_scores": scores.tolist(), "test_flagged": flagged.astype(int).tolist(),
        "TP": tp, "FP": fp, "FN": fn, "TN": tn,
        "TPR": tp / (tp + fn), "FPR": fp / (fp + tn),
        "precision": tp / (tp + fp) if tp + fp else None}


def main():
    torch.set_num_threads(1)
    torch.use_deterministic_algorithms(True)
    training = normal_points(160, 0.5)
    calibration = normal_points(40, 0.25)
    normal_test = normal_points(40, 0.75)
    abnormal_test = normal_points(20, 0.5)
    abnormal_test[:, 1] += np.where(np.arange(20) % 2 == 0, 0.8, -0.8)
    testing = np.vstack((normal_test, abnormal_test))
    labels = np.r_[np.zeros(40, dtype=int), np.ones(20, dtype=int)]
    center, scale = training.mean(axis=0), training.std(axis=0, ddof=0)
    train_scaled = (training - center) / scale
    calibration_scaled = (calibration - center) / scale
    test_scaled = (testing - center) / scale
    forest = IsolationForest(n_estimators=80, max_samples=128, random_state=SEED,
                             contamination="auto", n_jobs=1).fit(train_scaled)
    autoencoder, generator, training_record = fit_networks(
        torch.tensor(train_scaled, dtype=torch.float64))

    def score_models(values):
        x = torch.tensor(values, dtype=torch.float64)
        with torch.no_grad():
            reconstructed = autoencoder(x)
            _, z, z_reconstructed = generator(x)
            return {"isolation_forest": -forest.score_samples(values),
                "autoencoder_reconstruction_MSE": ((reconstructed-x)**2).mean(dim=1).numpy(),
                "EDE_latent_MSE": ((z-z_reconstructed)**2).mean(dim=1).numpy()}

    calibration_scores, test_scores = score_models(calibration_scaled), score_models(test_scaled)
    rows = [evaluate_scores(name, calibration_scores[name], test_scores[name], labels)
            for name in calibration_scores]
    for row in rows:
        assert np.isfinite(row["normal_calibration_scores"]).all()
        assert np.isfinite(row["test_scores"]).all()
        assert row["TP"] + row["FP"] + row["FN"] + row["TN"] == 60
    print(json.dumps({"scope": "Hand-set 2D manifold and displaced points; actual normal-only fits. "
        "EDE teaching architecture uses feature-MSE +50*contextual-L1 +latent-MSE, "
        "latent-MSE score like the retrieved official code, not exact paper equations/image replication.",
        "versions": {"python": platform.python_version(), "numpy": np.__version__,
            "sklearn": sklearn.__version__, "torch": torch.__version__},
        "device": "cpu", "dtype": "float64 for networks", "seed": SEED,
        "sizes": {"normal_training": 160, "normal_calibration": 40,
                  "normal_test": 40, "abnormal_test": 20},
        "train_only_mean": center.tolist(), "train_only_population_std": scale.tolist(),
        "training_record": training_record, "test_labels": labels.tolist(), "results": rows,
        "hypothetical_base_rate_example_not_model_measurement": {
            "assumed_TPR": 0.9, "assumed_FPR": 0.05, "assumed_prevalence": 0.01,
            "precision": 0.9*0.01/(0.9*0.01+0.05*0.99)},
        "limits": ["No abnormal labels or final test scores choose thresholds or training settings.",
            "The held-out normal quantile is empirical, not a probability or coverage guarantee.",
            "Deterministic synthetic sets are not iid external validation or a blind benchmark.",
            "Different raw score scales cannot be compared as probabilities.",
            "IF predict/auto offset is not used: the cutoff is calibrated on negative score_samples.",
            "EDE is a small fully connected concept model, not GANomaly image-benchmark reproduction.",
            "Paper training L2/inference L1 and retrieved software's MSE differ; no test-set minmax used.",
            "One seed and fixed budgets do not establish convergence or general model ranking.",
            "The artificial test anomaly fraction 20/60 is not a real operational base rate."]},
        ensure_ascii=False, indent=2))


if __name__ == "__main__":
    main()

假設異常比例 = 0.01
假設 TPR = 0.9，FPR = 0.05
Precision = (0.9 * 0.01) / (0.9 * 0.01 + 0.05 * 0.99)
          ≈ 0.153846

IMAGE AND DATA ATTRIBUTION

Figures 1–3 and Figure 5 are generated images from the official scikit-learn 1.9.1 documentation. Creator: The scikit-learn developers. Source pages: https://scikit-learn.org/stable/auto_examples/ensemble/plot_isolation_forest.html and https://scikit-learn.org/stable/auto_examples/model_selection/plot_precision_recall.html. The scikit-learn documentation FAQ states that its repository images and images generated in its documentation may be used under the BSD 3-Clause License. Figures 1–3 and Figure 5 are copied byte-for-byte from those pages.

Figure 4 is adapted from Figure 1 by Haoyang Jia and Wenfen Liu, “Anomaly detection in images with shared autoencoders,” published 2023-01-04, DOI: https://doi.org/10.3389/fnbot.2022.1046867, licensed CC BY 4.0: https://creativecommons.org/licenses/by/4.0/. The full source drawing was decoded from WebP and encoded as PNG; the mobile companion shows the full image followed by two overlapping source crops, resized for reading. This later simplified explanatory diagram represents the GANomaly method; the original method is by Samet Akcay, Amir Atapour-Abarghouei, and Toby P. Breckon, “GANomaly: Semi-Supervised Anomaly Detection via Adversarial Training” (2018): https://arxiv.org/abs/1805.06725. The original paper figure is not reproduced.

The precision-recall example uses scikit-learn’s `load_iris` data. Dataset credit: R. A. Fisher, “Iris” (1936), UCI Machine Learning Repository, DOI: https://doi.org/10.24432/C56C76. UCI lists the hosted Iris dataset under CC BY 4.0: https://creativecommons.org/licenses/by/4.0/. Scikit-learn documents that its v0.20+ packaged copy corrects two points according to Fisher and matches R, not the current UCI file. This corrected version and difference are disclosed; it is not claimed to be byte-identical to the current UCI download. The chart is a binary Iris/LinearSVC illustration, not anomaly-detection performance.

BSD 3-Clause License text (scikit-learn documentation images):

BSD 3-Clause License

Copyright (c) 2007-2026 The scikit-learn developers.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

* Neither the name of the copyright holder nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.