TimeRCD — 相対的文脈乖離による時系列異常検知基盤モデル【NeurIPS 2025】

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

def generate_normal_series(length=512):
    """正常な時系列を生成する（周期+トレンド+ノイズ）。"""
    t = np.arange(length)
    series = (np.sin(2 * np.pi * t / 50) +
              0.5 * np.sin(2 * np.pi * t / 120) +
              0.01 * t / length +
              np.random.randn(length) * 0.1)
    return series


def inject_anomalies(series, anomaly_type='spike', start=None, duration=10):
    """合成異常を注入する。"""
    s = series.copy()
    length = len(s)
    if start is None:
        start = length // 2

    labels = np.zeros(length)

    if anomaly_type == 'spike':
        spike_pos = start
        s[spike_pos] += 5 * np.std(series)
        labels[spike_pos] = 1

    elif anomaly_type == 'level_shift':
        s[start:start+duration] += 3 * np.std(series)
        labels[start:start+duration] = 1

    elif anomaly_type == 'trend_change':
        slope = 0.05
        for i in range(start, min(start+duration, length)):
            s[i] += slope * (i - start)
            labels[i] = 1

    elif anomaly_type == 'pattern_change':
        # 正弦波パターンを突然変更
        t_anom = np.arange(duration)
        s[start:start+duration] = 3 * np.sin(2 * np.pi * t_anom / 5)
        labels[start:start+duration] = 1

    return s, labels


def compute_rcd_scores(series, window_size=32, stride=1):
    """RCD（Relative Contextual Discrepancy）スコアを計算する。"""
    length = len(series)
    scores = np.zeros(length)

    for t in range(window_size, length - window_size, stride):
        left_window = series[t-window_size:t]
        right_window = series[t:t+window_size]

        # 各窓の統計的特徴を抽出（簡略化版）
        left_features = np.array([
            np.mean(left_window),
            np.std(left_window),
            np.mean(np.diff(left_window)),  # 平均変化率
            np.max(np.abs(np.diff(left_window))),  # 最大変化率
            np.polyfit(np.arange(window_size), left_window, 1)[0],  # トレンド
        ])

        right_features = np.array([
            np.mean(right_window),
            np.std(right_window),
            np.mean(np.diff(right_window)),
            np.max(np.abs(np.diff(right_window))),
            np.polyfit(np.arange(window_size), right_window, 1)[0],
        ])

        # ユークリッド距離
        scores[t] = np.sqrt(np.sum((left_features - right_features)**2))

    return scores


# 4種類の異常に対するRCDスコアを計算
anomaly_types = ['spike', 'level_shift', 'trend_change', 'pattern_change']
fig, axes = plt.subplots(4, 2, figsize=(16, 16))

for idx, atype in enumerate(anomaly_types):
    normal = generate_normal_series(512)
    anomalous, labels = inject_anomalies(normal, anomaly_type=atype,
                                          start=256, duration=30)
    rcd_scores = compute_rcd_scores(anomalous, window_size=32)

    # 左: 時系列
    axes[idx, 0].plot(normal, alpha=0.3, label='Normal', color='gray')
    axes[idx, 0].plot(anomalous, label='With anomaly', color='#00d4ff')
    # 異常区間をハイライト
    anomaly_mask = labels > 0
    if np.any(anomaly_mask):
        anom_start = np.where(anomaly_mask)[0][0]
        anom_end = np.where(anomaly_mask)[0][-1]
        axes[idx, 0].axvspan(anom_start, anom_end, alpha=0.2, color='red')
    axes[idx, 0].set_title(f'Time Series ({atype})')
    axes[idx, 0].set_ylabel('Value')
    axes[idx, 0].legend(fontsize=8)
    axes[idx, 0].grid(True, alpha=0.3)

    # 右: RCDスコア
    axes[idx, 1].plot(rcd_scores, color='#ffa726')
    if np.any(anomaly_mask):
        axes[idx, 1].axvspan(anom_start, anom_end, alpha=0.2, color='red')
    axes[idx, 1].set_title(f'RCD Score ({atype})')
    axes[idx, 1].set_ylabel('Score')
    axes[idx, 1].grid(True, alpha=0.3)

axes[-1, 0].set_xlabel('Time Step')
axes[-1, 1].set_xlabel('Time Step')

plt.tight_layout()
plt.savefig('timercd_anomaly_detection.png', dpi=150, bbox_inches='tight')
plt.show()

def generate_synthetic_training_data(n_samples=500, length=256):
    """合成異常付き学習データを生成する。"""
    X = []
    Y = []

    for _ in range(n_samples):
        series = generate_normal_series(length)

        if np.random.random() < 0.5:
            # 異常を注入
            atype = np.random.choice(['spike', 'level_shift',
                                      'trend_change', 'pattern_change'])
            start = np.random.randint(length // 4, 3 * length // 4)
            duration = np.random.randint(5, 40)
            series, labels = inject_anomalies(series, atype, start, duration)
        else:
            labels = np.zeros(length)

        X.append(series)
        Y.append(labels)

    return np.array(X), np.array(Y)


# 合成データ生成
X_train, Y_train = generate_synthetic_training_data(500, 256)

# RCDスコアと真のラベルの相関を評価
from sklearn.metrics import roc_auc_score

all_scores = []
all_labels = []

for i in range(len(X_train)):
    scores = compute_rcd_scores(X_train[i], window_size=24)
    all_scores.extend(scores[24:-24])  # 境界を除外
    all_labels.extend(Y_train[i][24:-24])

all_scores = np.array(all_scores)
all_labels = np.array(all_labels)

# 非ゼロスコアのみ
mask = all_scores > 0
if np.sum(all_labels[mask]) > 0:
    auc = roc_auc_score(all_labels[mask], all_scores[mask])
    print(f"RCD-based AUC-ROC: {auc:.4f}")

# 閾値による検出性能
thresholds = np.percentile(all_scores[mask], np.arange(80, 100, 1))
precisions = []
recalls = []

for thresh in thresholds:
    pred = (all_scores > thresh).astype(int)
    tp = np.sum((pred == 1) & (all_labels == 1))
    fp = np.sum((pred == 1) & (all_labels == 0))
    fn = np.sum((pred == 0) & (all_labels == 1))
    precision = tp / (tp + fp + 1e-10)
    recall = tp / (tp + fn + 1e-10)
    precisions.append(precision)
    recalls.append(recall)

plt.figure(figsize=(8, 5))
plt.plot(recalls, precisions, 'o-', color='#00d4ff', markersize=4)
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.title(f'RCD Anomaly Detection: Precision-Recall (AUC-ROC={auc:.3f})')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('timercd_pr_curve.png', dpi=150, bbox_inches='tight')
plt.show()