Anomaly Detection

import numpy as np

# 100,000 transactions, 100 of them fraudulent.
n, n_bad = 100_000, 100
prevalence = n_bad / n                     # 0.001

# A model that predicts 'normal' every single time:
accuracy = (n - n_bad) / n                 # 0.999
recall = 0.0                               # it catches nothing
print(f"prevalence {prevalence:.3%}  accuracy {accuracy:.3%}  recall {recall:.1%}")

# And a five-fold split leaves only about 20 positives per fold to evaluate on.
n_bad / 5

from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor

# Outlier detection: flag the odd rows in the data you have.
iso = IsolationForest(contamination=0.02, random_state=0)
flags = iso.fit_predict(X)                  # -1 anomaly, +1 normal

# LOF is transductive by default: fit_predict only, no predict.
lof = LocalOutlierFactor(n_neighbors=20)
flags = lof.fit_predict(X)
lof.negative_outlier_factor_                # the scores for the fitted rows

# Novelty detection: fit on clean data, then score new rows.
lof_nov = LocalOutlierFactor(n_neighbors=20, novelty=True).fit(X_clean)
lof_nov.predict(X_new)
lof_nov.decision_function(X_new)            # negative = anomalous

import numpy as np

x = np.array([50, 52, 48, 51, 49, 53, 47, 50, 52, 49,
              78, 51, 48, 50, 52, 49, 51, 47, 50, 22], dtype=float)

z = (x - x.mean()) / x.std(ddof=0)
flagged = np.abs(z) > 3
print(x.mean(), x.std(ddof=0))     # 49.95, 9.13
print(x[flagged], z[flagged])      # [78. 22.] [ 3.07 -3.06] — only just over

# Multivariate: the same idea with the covariance folded in.
from scipy.spatial.distance import mahalanobis

import numpy as np
from scipy.stats import median_abs_deviation

q1, q3 = np.percentile(x, [25, 75])
iqr = q3 - q1
lo, hi = q1 - 1.5 * iqr, q3 + 1.5 * iqr
tukey = (x < lo) | (x > hi)

mad = median_abs_deviation(x, scale=1.0)
mod_z = 0.6745 * (x - np.median(x)) / mad
robust = np.abs(mod_z) > 3.5

print(lo, hi, x[tukey])            # 44.5 56.5  [78. 22.]
print(np.round(mod_z[[10, 19]], 1))   # [ 12.8 -12.8]

import numpy as np
from sklearn.ensemble import IsolationForest

iso = IsolationForest(
    n_estimators=200,
    max_samples=256,          # each tree is grown on a small subsample, by design
    contamination=0.01,       # or "auto" to leave the threshold unset
    random_state=0,
).fit(X_train)

score = -iso.score_samples(X_test)     # higher = more anomalous
flags = iso.predict(X_test) == -1
threshold = -iso.offset_               # the cut contamination implied

# Rank by score rather than trusting contamination; review the top 200 by hand.
top = np.argsort(score)[::-1][:200]

from sklearn.svm import OneClassSVM, SGDOneClassSVM
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

oc = make_pipeline(
    StandardScaler(),
    OneClassSVM(kernel="rbf", nu=0.02, gamma="scale"),
).fit(X_clean)                 # fitted on data believed to be normal

oc.predict(X_new)              # -1 outside the boundary, +1 inside
oc.decision_function(X_new)    # signed distance: negative is outside

# Kernel OneClassSVM is quadratic-to-cubic in rows. Above ~20k rows:
SGDOneClassSVM(nu=0.02, random_state=0)

import numpy as np
from sklearn.neighbors import LocalOutlierFactor

lof = LocalOutlierFactor(n_neighbors=20, metric="euclidean")
flags = lof.fit_predict(X)                     # -1 anomaly, +1 normal
scores = -lof.negative_outlier_factor_         # ~1.0 normal, >1.5 suspicious

print(np.round(np.sort(scores)[::-1][:5], 2))

# For scoring unseen rows, refit in novelty mode on clean data.
nov = LocalOutlierFactor(n_neighbors=20, novelty=True).fit(X_clean)
nov.decision_function(X_new)

import numpy as np
from sklearn.decomposition import PCA

# The linear version: fit on clean rows, score by reconstruction error.
pca = PCA(n_components=10, random_state=0).fit(X_clean_scaled)
err = np.mean((X_new_scaled - pca.inverse_transform(pca.transform(X_new_scaled))) ** 2, axis=1)

# Threshold from the clean data's own error distribution, not from X_new.
clean_recon = pca.inverse_transform(pca.transform(X_clean_scaled))
clean_err = np.mean((X_clean_scaled - clean_recon) ** 2, axis=1)
threshold = float(np.quantile(clean_err, 0.995))
flags = err > threshold
print(f"threshold {threshold:.4f}, flagged {flags.mean():.2%}")

import numpy as np
from sklearn.metrics import (average_precision_score, classification_report,
                             precision_recall_curve, roc_auc_score)

score = -iso.score_samples(X_test)         # higher = more anomalous

roc_auc_score(y_test, score)               # flattering at low prevalence
average_precision_score(y_test, score)     # the number to quote instead

# precision@k: what the analyst reviewing 200 alerts a day actually experiences.
k = 200
top = np.argsort(score)[::-1][:k]
print(f"precision@{k} = {y_test[top].mean():.3f}")

print(classification_report(y_test, (score > threshold).astype(int), digits=3))

import numpy as np

# Capacity-driven: the team can review 50 alerts a day out of 40,000 rows.
score = -iso.score_samples(X_today)
budget = 50
threshold = float(np.sort(score)[-budget])
flags = score >= threshold

# Distribution-driven: a fixed quantile of a clean reference period.
threshold = float(np.quantile(reference_score, 0.999))

# Recheck it as the data drifts — the same cut flags a different fraction later.
for week, s in weekly_scores.items():
    print(week, f"{(s > threshold).mean():.3%} flagged")