Hands-on Lab: NYC Taxi Talep Anomalisinde 5 İstatistiksel Detektör Karşılaştırma

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import json
 
DATA_PATH = '../data/raw/numenta/NAB/data/realKnownCause/nyc_taxi.csv'
LABELS_PATH = '../data/raw/numenta/NAB/labels/combined_windows.json'
 
# Veriyi yükle
df = pd.read_csv(DATA_PATH, parse_dates=['timestamp'])
df = df.rename(columns={'value': 'passengers'})
 
print(f"Şekil: {df.shape}")
print(f"Zaman aralığı: {df['timestamp'].min()} → {df['timestamp'].max()}")
print(f"Frekans: {(df['timestamp'].iloc[1] - df['timestamp'].iloc[0])}")
df.head()

# Anomaly window'ları yükle
with open(LABELS_PATH) as f:
    all_labels = json.load(f)
 
nyc_windows = all_labels['realKnownCause/nyc_taxi.csv']
print(f"NYC Taxi anomaly window sayısı: {len(nyc_windows)}")
for w in nyc_windows:
    print(f"  {w[0]} → {w[1]}")
 
# Her satır için 'is_anomaly' etiketi üret
df['is_anomaly'] = 0
for start, end in nyc_windows:
    mask = (df['timestamp'] >= start) & (df['timestamp'] <= end)
    df.loc[mask, 'is_anomaly'] = 1
 
print(f"\nAnomaly satır sayısı: {df['is_anomaly'].sum()}")
print(f"Anomaly oranı: {df['is_anomaly'].mean()*100:.2f}%")

fig, ax = plt.subplots(figsize=(15, 5))
ax.plot(df['timestamp'], df['passengers'], color='steelblue',
        linewidth=0.5, alpha=0.7)
 
# Anomaly window'ları gölgele
for start, end in nyc_windows:
    ax.axvspan(start, end, color='crimson', alpha=0.2)
 
ax.set_title('NYC Taxi — Saatlik Yolcu Sayısı (kırmızı bölgeler anomaly)', fontsize=14)
ax.set_xlabel('Zaman')
ax.set_ylabel('Yolcu sayısı')
plt.tight_layout()
plt.savefig('reports/nyc_taxi_overview.png', dpi=120)
plt.show()
 
# İstatistikler
print(df['passengers'].describe())

from scipy.stats import genpareto
from statsmodels.stats.stattools import medcouple
 
def detector_zscore(x, threshold=3.0):
    """Klasik z-score detector."""
    z = (x - x.mean()) / (x.std(ddof=0) + 1e-9)
    return np.abs(z), np.abs(z) > threshold
 
def detector_modified_z(x, threshold=3.5):
    """Modified z-score (MAD-based)."""
    med = np.median(x)
    mad = np.median(np.abs(x - med))
    if mad == 0:
        return np.zeros_like(x), np.zeros_like(x, dtype=bool)
    m = 0.6745 * (x - med) / mad
    return np.abs(m), np.abs(m) > threshold
 
def detector_iqr(x, k=1.5):
    """IQR / Tukey's fences detector."""
    q1, q3 = np.percentile(x, [25, 75])
    iqr = q3 - q1
    lower, upper = q1 - k * iqr, q3 + k * iqr
    # Skor: sınırdan ne kadar uzakta?
    score = np.maximum(np.maximum(lower - x, 0), np.maximum(x - upper, 0))
    return score, (x < lower) | (x > upper)
 
def detector_adjusted_boxplot(x, k=1.5):
    """Hubert-Vandervieren adjusted boxplot."""
    q1, q3 = np.percentile(x, [25, 75])
    iqr = q3 - q1
    mc = medcouple(x)
    if mc >= 0:
        lower = q1 - k * np.exp(-4 * mc) * iqr
        upper = q3 + k * np.exp(3 * mc) * iqr
    else:
        lower = q1 - k * np.exp(-3 * mc) * iqr
        upper = q3 + k * np.exp(4 * mc) * iqr
    score = np.maximum(np.maximum(lower - x, 0), np.maximum(x - upper, 0))
    return score, (x < lower) | (x > upper)
 
def detector_pot(x, threshold_percentile=95, anomaly_q=0.001):
    """Peak Over Threshold with GPD."""
    u = np.percentile(x, threshold_percentile)
    excess = x[x > u] - u
    if len(excess) < 30:
        return np.zeros_like(x, dtype=float), np.zeros_like(x, dtype=bool)
    shape, _, scale = genpareto.fit(excess, floc=0)
 
    # Anomaly eşiği (VaR formulu)
    n, Nu = len(x), len(excess)
    if abs(shape) < 1e-6:
        var_th = u - scale * np.log((n / Nu) * anomaly_q)
    else:
        var_th = u + (scale / shape) * (((n / Nu) * anomaly_q) ** (-shape) - 1)
 
    # Skor: eşiği aşma miktarı
    score = np.maximum(x - u, 0)
    return score, x > var_th

from sklearn.metrics import (roc_auc_score, average_precision_score,
                              precision_score, recall_score, f1_score)
 
x = df['passengers'].values
y_true = df['is_anomaly'].values
 
detectors = {
    'Z-Score (k=3)':        lambda x: detector_zscore(x, 3.0),
    'Modified Z (k=3.5)':   lambda x: detector_modified_z(x, 3.5),
    'IQR (k=1.5)':          lambda x: detector_iqr(x, 1.5),
    'Adjusted Boxplot':     lambda x: detector_adjusted_boxplot(x, 1.5),
    'POT (q=0.001)':        lambda x: detector_pot(x, 95, 0.001),
}
 
results = []
for name, fn in detectors.items():
    scores, preds = fn(x)
    try:
        auc = roc_auc_score(y_true, scores)
        ap = average_precision_score(y_true, scores)
    except ValueError:
        auc, ap = np.nan, np.nan
    if preds.sum() > 0:
        prec = precision_score(y_true, preds, zero_division=0)
        rec = recall_score(y_true, preds, zero_division=0)
        f1 = f1_score(y_true, preds, zero_division=0)
    else:
        prec, rec, f1 = 0, 0, 0
    results.append({
        'Detector': name,
        'ROC-AUC': auc,
        'PR-AUC':  ap,
        'Precision': prec,
        'Recall':    rec,
        'F1':        f1,
        'Alarm rate': preds.mean(),
    })
 
results_df = pd.DataFrame(results)
print(results_df.round(3).to_string(index=False))

fig, axes = plt.subplots(len(detectors), 1, figsize=(16, 14), sharex=True)
 
for ax, (name, fn) in zip(axes, detectors.items()):
    scores, preds = fn(x)
    ax.plot(df['timestamp'], df['passengers'],
            color='steelblue', linewidth=0.3, alpha=0.6)
 
    # Tespit edilen anomaly'ler
    detected_mask = preds & ~y_true.astype(bool)  # tespit ama gerçek değil = FP
    correct_mask = preds & y_true.astype(bool)    # tespit + gerçek = TP
 
    ax.scatter(df.loc[detected_mask, 'timestamp'],
               df.loc[detected_mask, 'passengers'],
               c='orange', s=10, alpha=0.5, label='FP (tespit ama gerçek değil)')
    ax.scatter(df.loc[correct_mask, 'timestamp'],
               df.loc[correct_mask, 'passengers'],
               c='red', s=20, alpha=0.8, label='TP (doğru tespit)')
 
    # Gerçek anomaly window'ları gölgele
    for start, end in nyc_windows:
        ax.axvspan(start, end, color='crimson', alpha=0.08)
 
    ax.set_title(name, fontsize=11)
    ax.legend(loc='upper left', fontsize=8)
    ax.set_ylabel('Yolcu')
 
plt.tight_layout()
plt.savefig('reports/nyc_taxi_detectors.png', dpi=120)
plt.show()

def rolling_zscore(x, window=24):
    """24 saatlik rolling z-score (saatlik veri için günlük bağlam)."""
    s = pd.Series(x)
    mu = s.rolling(window, min_periods=window).mean()
    sigma = s.rolling(window, min_periods=window).std()
    return ((s - mu) / (sigma + 1e-9)).abs().fillna(0).values
 
scores_rolling = rolling_zscore(x, window=24)
preds_rolling = scores_rolling > 3.0
 
auc = roc_auc_score(y_true, scores_rolling)
ap = average_precision_score(y_true, scores_rolling)
prec = precision_score(y_true, preds_rolling, zero_division=0)
rec = recall_score(y_true, preds_rolling, zero_division=0)
f1 = f1_score(y_true, preds_rolling, zero_division=0)
 
print(f"Rolling Z-Score (24h)")
print(f"  ROC-AUC: {auc:.3f}")
print(f"  PR-AUC:  {ap:.3f}")
print(f"  P: {prec:.3f}, R: {rec:.3f}, F1: {f1:.3f}")
print(f"  Alarm rate: {preds_rolling.mean():.3f}")