導入
研究室でも現場でも、ログやセンサの時系列データは増え続けています。
- 実験装置の温度・圧力・流量
- 生体信号（心拍、SpO2、呼吸、加速度）
- 施設・都市の環境センサ
- アプリやサービスの監視指標（遅延、エラー率）
こうした時系列は「たまに壊れる」「たまにズレる」「たまに欠ける」一方で、全部を人間が監視するのは不可能です。
そこで重要になるのが **時系列異常検知（Anomaly Detection）**です。
本記事では、ダミーデータで「現場っぽい異常」を作り、
- 異常スコアを計算して
- 閾値でアラートを出し
- “どのセンサが怪しいか”の原因候補も出して
- 監視運用の落とし穴（誤検知・アラート疲れ・欠損・ドリフト）まで整理
する、Google Colab で動く最小例を作ります。
このデモは教育用です。特定分野（医療・材料・プロセス等）の真の分布や規格を再現するものではありません。
ただし「監視→検知→アラート→運用」という流れは多くの○○インフォマティクスで共通するため、方法論の理解に使えます。
TL;DR
- 異常検知は「モデル精度」より **運用（いつ・誰に・何を通知し、どう判断するか）**が勝負。
- 今回は3センサのダミー時系列に対し、学習データ（正常期間）だけで学習する Isolation Forest で異常スコアを出す。
- 閾値は「学習期間の上位1%」など 許容する誤検知率（false positive）から決めるのが実務的。
- アラートは点（1サンプル）ではなく、イベント（連続区間）にまとめ、クールダウンや最小継続時間でアラート疲れを減らす。
- “原因候補”は難しいが、まずは どのセンサがどれだけ基準から外れたかを出すだけでも運用が回りやすくなる。
1. まず「異常」の定義を決める（ここが一番大事）
異常検知の失敗は、たいてい「異常の定義」が曖昧なことから始まります。
本記事では分かりやすく、次の4タイプを異常として扱います。
- スパイク：瞬間的に跳ねる（ノイズより明らかに大きい）
- レベルシフト：平均が段差的に変わる（センサ校正ズレ等）
- ドリフト：じわじわ変わる（劣化・汚れ・季節変動の未吸収）
- 欠損：測定できていない（通信断・センサ停止）
現実でも、だいたいこのどれか（または混合）です。
2. Google Colabで動かす（コピペでOK）
以下のセルを 上から順に Colab に貼り付けて実行してください。
実行すると fig/ に画像が保存されます。
セル1：準備
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import IsolationForest
SEED = 42
rng = np.random.default_rng(SEED)
os.makedirs("fig", exist_ok=True)
セル2：ダミーの多変量時系列を生成（正常 + 異常 + 欠損）
def generate_multisensor_timeseries(n=2000, seed=42):
rng = np.random.default_rng(seed)
t = np.arange(n)
ts = pd.date_range("2025-01-01", periods=n, freq="min") # 1分刻み（仮）
# --- 正常時のベース（周期 + 小ノイズ） ---
s1 = 50 + 2.0*np.sin(2*np.pi*t/200) + rng.normal(0, 0.4, n) # Sensor1（例：温度）
s2 = 101 + 0.8*np.sin(2*np.pi*t/500 + 0.4) + rng.normal(0, 0.3, n) # Sensor2（例：圧力）
s3 = 0.20 + 0.05*np.sin(2*np.pi*t/50) + rng.normal(0, 0.02, n) # Sensor3（例：振動）
y = np.zeros(n, dtype=int) # 1なら異常（このデモでは「いずれかのセンサが異常」）
reason = np.array([""]*n, dtype=object)
def mark(idx, msg):
nonlocal y, reason
y[idx] = 1
reason[idx] = (reason[idx] + ("|" if reason[idx] else "") + msg)
# --- 異常1：スパイク（短いバースト） ---
for start, amp in [(750, 10.0), (1100, -12.0), (1700, 14.0)]:
idx = np.arange(start, min(start+4, n))
s1[idx] += amp
for i in idx:
mark(i, "spike(s1)")
# --- 異常2：レベルシフト（平均が変わる） ---
a, b = 900, 1040
s2[a:b] += 4.5
for i in range(a, b):
mark(i, "level_shift(s2)")
# --- 異常3：ドリフト（じわじわ変わる） ---
d0 = 1500
drift = 0.015*(t - d0)
drift[t < d0] = 0
s3 += drift * 0.3 # 小さめのドリフト（気づきにくい想定）
for i in range(d0, n):
mark(i, "drift(s3)")
# --- 異常4：欠損（通信断） ---
m0, m1 = 600, 650
s2[m0:m1] = np.nan
for i in range(m0, m1):
mark(i, "missing(s2)")
# ランダム欠損も少しだけ（現場っぽさ）
miss_rate = 0.01
for col in ["s1", "s2", "s3"]:
mask = rng.random(n) < miss_rate
if col == "s1":
s1[mask] = np.nan
for i in np.where(mask)[0]:
mark(i, "missing(s1)")
elif col == "s2":
s2[mask] = np.nan
for i in np.where(mask)[0]:
mark(i, "missing(s2)")
else:
s3[mask] = np.nan
for i in np.where(mask)[0]:
mark(i, "missing(s3)")
df = pd.DataFrame({
"time": ts,
"s1": s1,
"s2": s2,
"s3": s3,
"y_true": y,
"reason": reason,
})
return df
df = generate_multisensor_timeseries(n=2000, seed=SEED)
df.head()
セル3：図1（生データ + 真の異常区間）を描画して保存
def contiguous_intervals(mask_bool):
"""Trueが連続する区間[(start,end), ...]を返す（endは含む）"""
mask = np.asarray(mask_bool, dtype=bool)
n = len(mask)
intervals = []
i = 0
while i < n:
if not mask[i]:
i += 1
continue
s = i
while i < n and mask[i]:
i += 1
e = i - 1
intervals.append((s, e))
return intervals
def plot_raw_timeseries(df, savepath):
sensors = ["s1", "s2", "s3"]
intervals = contiguous_intervals(df["y_true"].values == 1)
fig, axes = plt.subplots(len(sensors), 1, figsize=(10, 7), sharex=True)
for ax, col in zip(axes, sensors):
ax.plot(df["time"], df[col])
ax.set_ylabel(col)
for s, e in intervals:
ax.axvspan(df["time"].iloc[s], df["time"].iloc[e], alpha=0.15)
ax.grid(True)
axes[-1].set_xlabel("time")
fig.suptitle("Raw time series (shaded = ground-truth anomalies in this toy demo)")
plt.tight_layout()
plt.savefig(savepath, dpi=160, bbox_inches="tight")
plt.show()
plot_raw_timeseries(df, "fig/fig1_raw_timeseries.png")
print("Saved:", "fig/fig1_raw_timeseries.png")
セル4：特徴量を作って学習（Isolation Forest）
ポイント：
- 実務ではセンサごとにスケールが違うので 標準化（StandardScaler）が効きます
- 欠損は運用上重要なので、欠損フラグも特徴量に入れます（「欠けたこと自体が異常」になり得るため）
SENSORS = ["s1", "s2", "s3"]
WINDOW = 25 # ローリング窓（短すぎるとノイズ、長すぎると鈍る）
TRAIN_END = 500 # ここまでは“正常運転で学習”したい想定（デモでは異常は主に後半に配置）
df_feat = df.copy()
# 欠損フラグ（0/1）
for col in SENSORS:
df_feat[f"{col}_isna"] = df_feat[col].isna().astype(int)
# 欠損値の埋め（簡易：補間→前方→後方）
for col in SENSORS:
df_feat[col] = df_feat[col].interpolate(limit_direction="both").ffill().bfill()
# 特徴量：値 / 差分 / ローリング平均 / ローリング標準偏差
for col in SENSORS:
df_feat[f"{col}_diff"] = df_feat[col].diff().fillna(0.0)
df_feat[f"{col}_mean"] = df_feat[col].rolling(WINDOW, min_periods=1).mean()
df_feat[f"{col}_std"] = df_feat[col].rolling(WINDOW, min_periods=1).std(ddof=0).fillna(0.0)
feature_cols = []
for col in SENSORS:
feature_cols += [col, f"{col}_diff", f"{col}_mean", f"{col}_std", f"{col}_isna"]
X = df_feat[feature_cols].values
y_true = df_feat["y_true"].values
# 標準化
scaler = StandardScaler()
X_train = X[:TRAIN_END]
X_train_s = scaler.fit_transform(X_train)
X_all_s = scaler.transform(X)
# Isolation Forest（教師なし）
iso = IsolationForest(
n_estimators=300,
random_state=SEED,
contamination="auto",
)
iso.fit(X_train_s)
# 異常スコア（大きいほど異常にする）
score_train = -iso.decision_function(X_train_s)
score_all = -iso.decision_function(X_all_s)
print("score_train range:", float(score_train.min()), float(score_train.max()))
print("score_all range:", float(score_all.min()), float(score_all.max()))
セル5：アラート化（閾値・最小継続・クールダウン）
ここが「監視と運用」のコアです。
- 閾値：学習期間のスコア上位
p%を超えたらアラート
→ “許容する誤検知率”から決められる - 最小継続：短いノイズで鳴らさない
- クールダウン：同一事象で鳴り続けない（アラート疲れ対策）
# 閾値（例：学習期間で上位1%を超えたら異常）
Q = 0.99
threshold = float(np.quantile(score_train, Q))
alert_raw = score_all > threshold
def extract_events(alert_flags, score, times, min_len=3, cooldown=20):
"""アラートの連続区間をイベントにまとめる（簡易）"""
n = len(alert_flags)
events = []
i = 0
eid = 0
while i < n:
if not alert_flags[i]:
i += 1
continue
s = i
while i < n and alert_flags[i]:
i += 1
e = i - 1
length = e - s + 1
if length >= min_len:
seg = score[s:e+1]
peak_rel = int(np.argmax(seg))
peak = s + peak_rel
events.append({
"event_id": eid,
"start": times[s],
"end": times[e],
"length": length,
"peak_time": times[peak],
"peak_score": float(score[peak]),
"peak_idx": peak,
})
eid += 1
i = e + 1 + cooldown
else:
i = e + 1
return events
events = extract_events(alert_raw, score_all, df_feat["time"].values, min_len=3, cooldown=20)
print("threshold =", threshold)
print("n_alert_points =", int(alert_raw.sum()))
print("n_events =", len(events))
pd.DataFrame(events).head()
セル6：図2（異常スコア + 閾値 + 真の異常区間）を描画
def plot_score(df, score, threshold, alert_flags, savepath):
intervals = contiguous_intervals(df["y_true"].values == 1)
plt.figure(figsize=(10,4.5))
plt.plot(df["time"], score, label="anomaly score")
plt.axhline(threshold, linestyle="--", label=f"threshold (train Q={Q})")
# 真の異常区間（このデモ用）
for s, e in intervals:
plt.axvspan(df["time"].iloc[s], df["time"].iloc[e], alpha=0.12)
# アラート点（散布）
idx = np.where(alert_flags)[0]
plt.scatter(df["time"].iloc[idx], score[idx], s=12)
plt.title("Anomaly score & alerts (shaded = ground-truth anomalies in this toy demo)")
plt.xlabel("time")
plt.ylabel("score (higher = more anomalous)")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(savepath, dpi=160, bbox_inches="tight")
plt.show()
plot_score(df_feat, score_all, threshold, alert_raw, "fig/fig2_score_and_alerts.png")
print("Saved:", "fig/fig2_score_and_alerts.png")
セル7：原因候補（どのセンサがどれだけ外れた？）を出す + 図3
ここでは、厳密な因果はやりません。運用でまず欲しいのは「どのセンサが目立ってズレたか」です。
- 各センサについて「ローリング平均との差」をローリング標準偏差で割った zスコア（簡易）を作る
- イベントのピーク時点で、|z| が大きいセンサを原因候補として出す
eps = 1e-9
z = pd.DataFrame(index=df_feat.index)
for col in SENSORS:
resid = df_feat[col] - df_feat[f"{col}_mean"]
z[col] = resid / (df_feat[f"{col}_std"] + eps)
def top_sensors_at(idx, z_df, k=2):
s = z_df.loc[idx, SENSORS].abs().sort_values(ascending=False)
return list(s.index[:k]), list(s.values[:k])
# イベント表を作る（上位2センサを付与）
rows = []
for ev in events:
idx = ev["peak_idx"]
tops, vals = top_sensors_at(idx, z)
rows.append({
"event_id": ev["event_id"],
"start": ev["start"],
"end": ev["end"],
"length": ev["length"],
"peak_time": ev["peak_time"],
"peak_score": ev["peak_score"],
"top_sensors": ", ".join([f"{s}(|z|={v:.1f})" for s, v in zip(tops, vals)]),
"reason_gt": df_feat.loc[idx, "reason"], # デモなので真の理由も見える
})
event_df = pd.DataFrame(rows).sort_values("peak_score", ascending=False)
event_df.head(10)
# 最も大きいイベントで「原因候補」を棒グラフにする
if len(event_df) > 0:
top_ev = event_df.iloc[0]
peak_idx = int(events[int(top_ev["event_id"])]["peak_idx"])
contrib = z.loc[peak_idx, SENSORS].abs()
plt.figure(figsize=(6,4))
plt.bar(contrib.index, contrib.values)
plt.title("Root-cause hint (simple): |z-score| at event peak")
plt.ylabel("|z| (deviation from rolling baseline)")
plt.tight_layout()
plt.savefig("fig/fig3_rootcause_hint.png", dpi=160, bbox_inches="tight")
plt.show()
print("Saved:", "fig/fig3_rootcause_hint.png")
else:
print("No events detected. Try lowering Q (e.g., 0.98) or reducing min_len.")
セル8：おまけ（評価：このデモでは真のラベルがあるので計算できる）
現実ではラベルが無いことが多いですが、デモでは「どれくらい当たったか」を確認できます。
# 点（時刻）レベルの簡易評価
tp = int(np.sum((alert_raw == 1) & (y_true == 1)))
fp = int(np.sum((alert_raw == 1) & (y_true == 0)))
fn = int(np.sum((alert_raw == 0) & (y_true == 1)))
precision = tp / (tp + fp + 1e-12)
recall = tp / (tp + fn + 1e-12)
f1 = 2 * precision * recall / (precision + recall + 1e-12)
print("Point-level precision =", precision)
print("Point-level recall =", recall)
print("Point-level F1 =", f1)
print("TP, FP, FN =", tp, fp, fn)
3. 本記事掲載コードを実行して得られる図
図1：生データ（真の異常区間を影で表示）
図2：異常スコアとアラート（閾値超え）
図3：原因候補（どのセンサが外れたか）
4. 図の読み方（特に図2と図3）
図2：スコアが閾値を超えたら「アラート」
- 折れ線：異常スコア（高いほど異常っぽい）
- 破線：閾値（学習期間の上位1% など）
- 点：閾値を超えた時刻（alert）
- 影：真の異常区間（このデモでだけ見える答え）
実務では「影（真の異常）」は見えないので、閾値・アラート疲れ・原因候補の設計が重要になります。
図3：原因“候補”としての |z|（簡易）
図3は「イベントのピーク時刻に、どのセンサがローリング平均との差から大きく外れたか」を棒グラフ化しています。
- これは因果推論ではありません
- でも運用上は、
“どのセンサから見ればよいか”の第一候補として役に立ちます
今回の閾値は「学習期間の異常スコア上位1%（Q=0.99）」に設定しており、誤検知を抑える“保守的”な設定です。
そのため、点（時刻）単位の評価では Precision は高くなりやすい一方、長い異常区間（例：ドリフト）では Recall が低く見えることがあります。
監視運用では「点」よりも「イベントとして一度でも鳴ったか」「どれくらい早く鳴ったか（検知遅延）」が重要になることが多いです。
5. 監視・運用でハマりやすいポイント（ここが本題）
5.1 閾値は「精度」より「許容誤検知率」で決めると回りやすい
現場では「誤検知ゼロ」はほぼ不可能です。
むしろ 1日あたり何件までなら対応できるか が重要です。
このデモのように
- 学習期間の上位
1%を閾値にする - あるいは「1時間に平均1回まで」など運用都合で決める
という考え方が実務的です。
5.2 アラートは“点”ではなく“イベント”にまとめる
1分おきに鳴ったら誰も見ません。
- 最小継続時間（min_len）
- クールダウン（cooldown）
で アラート疲れを減らして、見るべきイベントに絞ります。
5.3 欠損は「異常」でもあり「前処理」でもある
欠損が起きたら、まず運用としては
- センサ停止なのか
- 通信断なのか
- 単なる欠測なのか
を切り分けたいことが多いです。
そのため本記事では、欠損フラグを特徴量に入れました。
5.4 ドリフトは最難関（じわじわ壊れる）
スパイクは比較的検知しやすいですが、ドリフトは難しいです。
- ベースラインが季節変動していると誤検知しやすい
- モデルも古くなる
対策は「再学習」だけではなく、
- モデルの監視（ドリフト検知）
- 参照期間の更新
- 目的の再定義（何が異常なのか）
が必要になります。
6. 研究室運営・教育に使うなら（発展課題）
- 閾値Qを変えて（0.98/0.995など）誤検知と見逃しのトレードオフを観察する
- 欠損を「異常」ではなく「データ品質問題」として別系統で通知する設計に変える
- イベント検出を「k out of m（m点中k点が異常なら通知）」に変えて頑健化する
- 予測モデル（1-step予測）を作り、予測誤差で異常検知する（より時系列らしい）
- 多変量の“原因推定”を、寄与度（簡易SHAPなど）やグラフ構造で改善する
まとめ
本記事では、ダミーデータで「監視の現場っぽい異常」を作り、
- 多変量時系列から特徴量を作る
- Isolation Forestで異常スコアを出す
- 閾値でアラート化し、イベントにまとめる
- 原因候補として「どのセンサが外れたか」を提示する
という 監視→検知→運用までの最小構成を、Google Colabで再現しました。
時系列異常検知は、モデル選び以上に
- 異常定義
- 閾値設計
- アラート疲れ対策
- 欠損・ドリフトへの対応
- 原因候補の提示
といった「運用設計」が成果を左右します。
まずはこの最小例を動かし、あなたの分野のセンサ・ログに置き換えてみてください。