センサー単体では見逃す設備異常を捉える：多変量異常検知の実務（No.061〜No.070）

製造設備では、温度や圧力がそれぞれ管理範囲内でも、複数センサーの関係が通常運転から崩れることがあります。本記事では、架空の連続生産設備を題材に、相関の確認、主成分分析（PCA）、再構成誤差、マハラノビス距離、原因候補の整理までを一つの流れで扱います。

対象は「異常検知100本ノック」の No.061〜No.070（第7章：多変量異常検知） です。目的はモデルを動かすことだけではなく、保全・品質・生産の担当者が、点検優先度や追加計測の判断に使える形へ落とし込むことです。

[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。

はじめに：この記事で扱う製造業の実務課題

連続生産設備「Mixer-A」では、温度・圧力・振動・モーター電流・冷却水流量を1分ごとに監視しています。各センサーの上下限監視だけでは、負荷に対して冷却が不足している、圧力と流量の釣り合いが崩れている、といった関係性の異常を見逃すことがあります。

今回の業務上の問いは、次の3点です。

1. 通常運転のセンサー間関係を、どのように数値化するか
2. 複数の検知法から、点検対象をどのように絞るか
3. 異常スコアを、原因調査の初動へどうつなぐか

現場でよくある状況

• センサーごとの管理限界は守られているのに、品質ばらつきが増える
• 設備負荷や品種の違いで正常範囲が動き、固定しきい値では誤報が増える
• センサー数を増やした結果、相関の確認や原因説明が難しくなる
• 異常スコアは出せても、保全担当者が「どこを確認するか」を決められない

多変量異常検知は、単変量監視を置き換えるものではありません。安全上限・下限は従来どおり監視し、その内側で生じる関係崩れを補完的に捉える位置づけが現実的です。

なぜこの問題は判断が難しいのか

センサーが $p$ 個ある観測を $\mathbf{x}\in\mathbb{R}^p$ とすると、正常・異常の境界は各軸に平行な箱ではなく、変数間の相関を反映した斜めの領域になります。また、稼働条件、品種、季節、センサー交換によって正常状態そのものも変化します。

そのため、実務では次を分けて考える必要があります。

• 統計的な逸脱：過去の通常データから離れているか
• 業務上の異常：安全・品質・停止リスクにつながるか
• 対応可能性：現場が確認できる原因候補に変換できるか

モデルのスコアは診断結果ではなく、点検順序を決めるための証拠の一つです。

今回扱うノックの全体像

No.	テーマ	実務上の成果物
061	多変量異常	単変量監視との差の理解
062	相関行列	通常時の関係性一覧
063	散布図行列	非線形・群・外れ方の目視確認
064	PCA	少数の運転状態軸への要約
065	PCA再構成誤差	通常構造からの外れ度
066	主成分空間	点検候補の可視化
067	マハラノビス距離	相関を考慮した距離
068	距離による判定	管理可能なアラート候補
069	原因候補分析	寄与の大きいセンサー候補
070	高次元の注意点	導入前チェックリスト

Python環境の準備

NumPy、pandas、matplotlib、scikit-learnを使用します。外部データは読み込まず、乱数シードを固定して再現可能にします。グラフ内の表記は実行環境による日本語フォント差を避けるため英語とします。

import sys
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import sklearn
from pandas.plotting import scatter_matrix
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

SEED = 202507
rng = np.random.default_rng(SEED)
pd.set_option("display.max_columns", 20)
pd.set_option("display.float_format", lambda x: f"{x:,.3f}")

print(f"Python: {sys.version.split()[0]}")
print(f"NumPy: {np.__version__}, pandas: {pd.__version__}")
print(f"matplotlib: {matplotlib.__version__}, scikit-learn: {sklearn.__version__}")

Python: 3.13.1
NumPy: 2.5.1, pandas: 3.0.3
matplotlib: 3.11.0, scikit-learn: 1.9.0

架空データの作成

480分の運転データを作ります。先頭360分を設備が安定していた基準期間、後半120分を監視期間とします。監視期間の一部に、次の関係崩れを注入します。

• 冷却系異常：負荷に対して温度が高く、冷却水流量が低い
• 駆動系異常：負荷に対して振動と電流が高い
• 圧送系異常：圧力と流量の通常関係が崩れる

true_cause は検証用の架空ラベルです。学習やしきい値決定には使用しません。

n = 480
time = pd.date_range("2025-01-15 08:00", periods=n, freq="min")
load = np.clip(rng.normal(0, 1, n), -2.4, 2.4)
ambient = 0.7 * np.sin(np.linspace(0, 3 * np.pi, n)) + rng.normal(0, 0.15, n)

df = pd.DataFrame({
    "timestamp": time,
    "temperature_C": 68.0 + 2.2 * load + 0.8 * ambient + rng.normal(0, 0.45, n),
    "pressure_MPa": 1.80 + 0.18 * load + rng.normal(0, 0.035, n),
    "vibration_mm_s": 2.40 + 0.28 * load + rng.normal(0, 0.10, n),
    "motor_current_A": 42.0 + 4.8 * load + rng.normal(0, 0.75, n),
    "coolant_flow_L_min": 31.0 + 2.0 * load - 0.5 * ambient + rng.normal(0, 0.55, n),
})
df["true_cause"] = "normal"

anomaly_map = {
    "cooling_relation": [386, 387, 388, 389],
    "drive_relation": [421, 422, 423, 424],
    "pressure_relation": [458, 459, 460, 461],
}
for cause, idx in anomaly_map.items():
    df.loc[idx, "true_cause"] = cause
df.loc[anomaly_map["cooling_relation"], ["temperature_C", "coolant_flow_L_min"]] += [3.2, -3.0]
df.loc[anomaly_map["drive_relation"], ["vibration_mm_s", "motor_current_A"]] += [0.75, 5.5]
df.loc[anomaly_map["pressure_relation"], ["pressure_MPa", "coolant_flow_L_min"]] += [0.22, -2.6]

features = ["temperature_C", "pressure_MPa", "vibration_mm_s", "motor_current_A", "coolant_flow_L_min"]
train_mask = np.arange(n) < 360
print(f"データ件数: {len(df):,}、基準期間: {train_mask.sum():,}、監視期間: {(~train_mask).sum():,}")
df.head()

データ件数: 480、基準期間: 360、監視期間: 120

	timestamp	temperature_C	pressure_MPa	vibration_mm_s	motor_current_A	coolant_flow_L_min	true_cause
0	2025-01-15 08:00:00	63.187	1.421	1.781	29.779	27.232	normal
1	2025-01-15 08:01:00	65.697	1.655	2.311	38.384	30.303	normal
2	2025-01-15 08:02:00	69.995	1.975	2.568	45.791	32.274	normal
3	2025-01-15 08:03:00	68.658	1.932	2.590	42.751	32.041	normal
4	2025-01-15 08:04:00	66.395	1.642	2.153	37.482	30.070	normal

No.061：多変量異常とは何かを理解する

実務での意味

上下限内というだけでは、設備の整合性までは保証できません。たとえば温度72℃と冷却水流量29 L/minは単独では許容範囲でも、高負荷時の通常関係から見れば不自然かもしれません。

分析・モデル化の考え方

単変量ルールを「温度 63〜75℃、圧力 1.4〜2.2 MPa、振動 1.6〜3.4 mm/s、電流 30〜55 A、流量 25〜38 L/min」と仮定し、既知の関係異常がすり抜けるか確認します。

Pythonで確認する

limits = {
    "temperature_C": (63, 75), "pressure_MPa": (1.4, 2.2),
    "vibration_mm_s": (1.6, 3.4), "motor_current_A": (30, 55),
    "coolant_flow_L_min": (25, 38),
}
single_alarm = pd.Series(False, index=df.index)
for col, (low, high) in limits.items():
    single_alarm |= ~df[col].between(low, high)

comparison = pd.crosstab(df["true_cause"], single_alarm, rownames=["true_cause"], colnames=["single_limit_alarm"])
comparison

single_limit_alarm	False	True
true_cause
cooling_relation	4	0
drive_relation	4	0
normal	453	15
pressure_relation	2	2

結果の読み取り

単変量ルールで検知できる点があっても、関係異常の全件を安定して拾えるとは限りません。逆に通常運転の端を異常とする可能性もあります。したがって、安全限界の監視は残しつつ、相関構造からの逸脱を別のスコアとして併用します。

No.062：相関行列を確認する

実務での意味

相関行列は、通常時にどのセンサーが同じ運転負荷に反応しているかを一覧化します。強い相関は冗長性の手掛かりになりますが、因果関係を証明するものではありません。

分析・モデル化の考え方

基準期間だけで Pearson の相関係数 $r$ を求めます。監視期間を混ぜると、異常そのものが「通常の関係」を変えてしまうためです。

Pythonで確認する

corr = df.loc[train_mask, features].corr()
corr.round(2)

	temperature_C	pressure_MPa	vibration_mm_s	motor_current_A	coolant_flow_L_min
temperature_C	1.000	0.940	0.910	0.950	0.900
pressure_MPa	0.940	1.000	0.930	0.970	0.930
vibration_mm_s	0.910	0.930	1.000	0.930	0.910
motor_current_A	0.950	0.970	0.930	1.000	0.940
coolant_flow_L_min	0.900	0.930	0.910	0.940	1.000

fig, ax = plt.subplots(figsize=(7, 5))
im = ax.imshow(corr, vmin=-1, vmax=1, cmap="coolwarm")
ax.set_xticks(range(len(features)), [c.replace("_", "\n") for c in features], rotation=35, ha="right")
ax.set_yticks(range(len(features)), [c.replace("_", " ") for c in features])
for i in range(len(features)):
    for j in range(len(features)):
        ax.text(j, i, f"{corr.iloc[i, j]:.2f}", ha="center", va="center", fontsize=8)
ax.set_title("Correlation matrix during baseline operation")
ax.set_xlabel("Sensor variable")
ax.set_ylabel("Sensor variable")
ax.grid(False)
fig.colorbar(im, ax=ax, label="Correlation coefficient")
plt.tight_layout()
plt.show()

結果の読み取り

電流・温度・圧力などに正の相関が見られるのは、共通の運転負荷に反応するよう生成したためです。ここで重要なのは係数の大きさそのものより、この基準関係が今後も維持されるかです。品種や運転モードが複数ある場合は、モード別に相関を確認します。

No.063：散布図行列で変数間の関係を見る

実務での意味

相関係数だけでは、曲線関係、複数の運転群、少数の外れ点を見落とします。散布図行列は、モデル構築前の現場ヒアリング対象を見つける診断図です。

分析・モデル化の考え方

表示を読みやすくするため、代表4変数について基準期間と監視期間の異常点を重ねます。色は真の原因ではなく、基準・監視正常・検証用異常の区分です。

Pythonで確認する

plot_cols = ["temperature_C", "pressure_MPa", "motor_current_A", "coolant_flow_L_min"]
view = df.loc[::3, plot_cols].copy()
axes = scatter_matrix(view, figsize=(9, 9), diagonal="hist", alpha=0.45, color="#4C78A8")
abnormal = df[df["true_cause"] != "normal"]
for i, y in enumerate(plot_cols):
    for j, x in enumerate(plot_cols):
        ax = axes[i, j]
        if i != j:
            ax.scatter(abnormal[x], abnormal[y], color="#D62728", s=22, marker="x")
        ax.grid(True, alpha=0.25)
plt.suptitle("Scatter matrix: baseline pattern and injected relation anomalies", y=1.01)
plt.xlabel("Sensor value / red x: injected anomaly")
plt.ylabel("Sensor value")
plt.tight_layout()
plt.show()

結果の読み取り

赤い点は、個々の軸では極端でなくても、通常の点群が作る細長い帯から外れます。群が複数に分かれる場合は、異常と決めつけず、品種・段取り・停止復帰などの運転モードを確認する必要があります。

No.064：PCAで次元削減する

実務での意味

PCAは、相関する多数のセンサーを少数の「運転状態軸」に要約します。監視画面の軸数を減らし、通常運転の主要な変動と、それ以外の小さな関係崩れを分けて考えられます。

分析・モデル化の考え方

単位の異なるセンサーを標準化し、基準期間だけでPCAを学習します。標準化後の行列 $\mathbf{Z}$ を、主成分負荷量 $\mathbf{P}$ で

$$\mathbf{T}=\mathbf{Z}\mathbf{P}$$

へ写像します。今回は可視化と再構成のため2主成分を採用しますが、採用数は累積寄与率だけでなく、見逃し・誤報コストも含めて決めます。

Pythonで確認する

scaler = StandardScaler()
X_train = scaler.fit_transform(df.loc[train_mask, features])
X_all = scaler.transform(df[features])

pca_full = PCA().fit(X_train)
explained = pd.DataFrame({
    "component": [f"PC{i}" for i in range(1, len(features) + 1)],
    "explained_ratio": pca_full.explained_variance_ratio_,
    "cumulative_ratio": np.cumsum(pca_full.explained_variance_ratio_),
})
explained

	component	explained_ratio	cumulative_ratio
0	PC1	0.945	0.945
1	PC2	0.020	0.965
2	PC3	0.019	0.984
3	PC4	0.011	0.994
4	PC5	0.006	1.000

fig, ax = plt.subplots(figsize=(7, 4))
ax.bar(explained["component"], explained["explained_ratio"], color="#4C78A8", label="Individual")
ax.plot(explained["component"], explained["cumulative_ratio"], color="#F58518", marker="o", label="Cumulative")
ax.set_title("PCA explained variance during baseline operation")
ax.set_xlabel("Principal component")
ax.set_ylabel("Explained variance ratio")
ax.set_ylim(0, 1.05)
ax.grid(True, axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

結果の読み取り

第1主成分が共通負荷を大きく表し、後続成分がセンサー間の細かな差を表します。寄与率の低い成分にも異常情報が残ることがあるため、「高い累積寄与率になるまで残す」だけではなく、再構成誤差の検証が必要です。

No.065：PCAの再構成誤差で異常検知する

実務での意味

通常運転を2主成分で再現しにくい観測は、主要な運転パターンから外れた点検候補です。設備負荷が高いだけの点と、センサー関係が崩れた点を分けやすくなります。

分析・モデル化の考え方

標準化観測 $\mathbf{z}_i$ と再構成値 $\hat{\mathbf{z}}_i$ の二乗誤差を

$$Q_i=\lVert\mathbf{z}_i-\hat{\mathbf{z}}_i\rVert_2^2$$

とします。しきい値は基準期間の99パーセンタイルとし、通常時に想定されるアラート率を明示します。

Pythonで確認する

pca = PCA(n_components=2).fit(X_train)
scores = pca.transform(X_all)
X_reconstructed = pca.inverse_transform(scores)
reconstruction_error = np.sum((X_all - X_reconstructed) ** 2, axis=1)
pca_threshold = np.quantile(reconstruction_error[train_mask], 0.99)
df["pca_error"] = reconstruction_error
df["pca_alarm"] = df["pca_error"] > pca_threshold

pd.DataFrame({
    "baseline_99pct_threshold": [pca_threshold],
    "monitoring_alarms": [int(df.loc[~train_mask, "pca_alarm"].sum())],
    "monitoring_alarm_rate": [df.loc[~train_mask, "pca_alarm"].mean()],
})

	baseline_99pct_threshold	monitoring_alarms	monitoring_alarm_rate
0	0.708	14	0.117

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(df["timestamp"], df["pca_error"], color="#4C78A8", linewidth=1, label="PCA reconstruction error")
ax.axhline(pca_threshold, color="#D62728", linestyle="--", label="Baseline 99th percentile")
ax.axvline(df.loc[360, "timestamp"], color="gray", linestyle=":", label="Monitoring start")
ax.scatter(df.loc[df["pca_alarm"], "timestamp"], df.loc[df["pca_alarm"], "pca_error"], color="#D62728", s=20)
ax.set_title("PCA reconstruction error over time")
ax.set_xlabel("Timestamp")
ax.set_ylabel("Squared reconstruction error")
ax.grid(True, alpha=0.3)
ax.legend(loc="upper left")
plt.tight_layout()
plt.show()

結果の読み取り

監視期間でしきい値を超えた時刻は、通常構造では説明しにくい候補です。99パーセンタイルは万能値ではありません。1分ごとの誤報率が小さくても、24時間運転では件数が積み上がるため、連続回数や設備停止コストを加味して運用しきい値を調整します。

No.066：主成分空間で異常点を可視化する

実務での意味

主成分空間は、アラートが「負荷の端にあるのか」「通常群から別方向へ外れたのか」を現場と共有する図になります。同じ方向に固まるアラートは、共通原因や連続事象の可能性があります。

分析・モデル化の考え方

横軸をPC1、縦軸をPC2とし、PCA再構成誤差のアラートを重ねます。ただし、2次元図で重なる点でも、捨てた成分方向には離れている場合があります。

Pythonで確認する

fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(scores[train_mask, 0], scores[train_mask, 1], s=15, alpha=0.35, color="#4C78A8", label="Baseline")
monitor_normal = (~train_mask) & (~df["pca_alarm"].to_numpy())
ax.scatter(scores[monitor_normal, 0], scores[monitor_normal, 1], s=18, alpha=0.55, color="#54A24B", label="Monitoring / no alarm")
alarm = df["pca_alarm"].to_numpy()
ax.scatter(scores[alarm, 0], scores[alarm, 1], s=42, marker="x", color="#D62728", label="PCA alarm")
ax.set_title("Operating states in principal-component space")
ax.set_xlabel("PC1 score")
ax.set_ylabel("PC2 score")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

結果の読み取り

点の位置は運転状態の類似性、色は判定結果を表します。主成分空間の端にあるだけで故障とは限りません。高負荷・立上げなどの運転履歴と照合し、同じ領域のアラートが繰り返すかを確認します。

No.067：マハラノビス距離を計算する

実務での意味

マハラノビス距離は、変数のばらつきと相関を考慮して「通常中心からどれだけ離れたか」を測ります。相関方向に沿った変化は比較的小さく、関係を横切る変化は大きく評価されます。

分析・モデル化の考え方

基準期間の平均 $\boldsymbol{\mu}$ と共分散行列 $\mathbf{S}$ に対し、二乗マハラノビス距離を

$$D_i^2=(\mathbf{z}_i-\boldsymbol{\mu})^\top\mathbf{S}^{-1}(\mathbf{z}_i-\boldsymbol{\mu})$$

とします。数値安定性のため、逆行列には擬似逆行列を使います。

Pythonで確認する

mu = X_train.mean(axis=0)
cov = np.cov(X_train, rowvar=False)
precision = np.linalg.pinv(cov)
centered = X_all - mu
mahalanobis_sq = np.einsum("ij,jk,ik->i", centered, precision, centered)
df["mahalanobis_sq"] = mahalanobis_sq

df.loc[~train_mask, ["timestamp", "true_cause", "mahalanobis_sq"]].nlargest(10, "mahalanobis_sq")

	timestamp	true_cause	mahalanobis_sq
424	2025-01-15 15:04:00	drive_relation	122.452
422	2025-01-15 15:02:00	drive_relation	108.729
421	2025-01-15 15:01:00	drive_relation	107.022
458	2025-01-15 15:38:00	pressure_relation	72.714
461	2025-01-15 15:41:00	pressure_relation	66.094
388	2025-01-15 14:28:00	cooling_relation	63.051
460	2025-01-15 15:40:00	pressure_relation	62.694
423	2025-01-15 15:03:00	drive_relation	60.361
459	2025-01-15 15:39:00	pressure_relation	57.027
387	2025-01-15 14:27:00	cooling_relation	53.034

結果の読み取り

上位一覧は、監視期間の点検優先順位として利用できます。ただし、距離が大きい理由は故障だけでなく、未学習の正常な品種・運転条件でも起こります。距離の順位と運転コンテキストをセットで保存します。

No.068：マハラノビス距離で異常検知する

実務での意味

距離をアラートへ変換すると、日次の点検件数や対応工数を計画できます。しきい値は統計的有意性だけでなく、現場が処理できる件数と見逃しコストで決めます。

分析・モデル化の考え方

分布仮定に依存しすぎないよう、ここでも基準期間の経験的99パーセンタイルを採用します。PCA再構成誤差との一致・不一致も確認し、異なる種類の逸脱を区別します。

Pythonで確認する

md_threshold = np.quantile(df.loc[train_mask, "mahalanobis_sq"], 0.99)
df["md_alarm"] = df["mahalanobis_sq"] > md_threshold

alarm_compare = pd.crosstab(
    df.loc[~train_mask, "pca_alarm"],
    df.loc[~train_mask, "md_alarm"],
    rownames=["PCA alarm"], colnames=["Mahalanobis alarm"],
)
print(f"Mahalanobis threshold (baseline 99th percentile): {md_threshold:.3f}")
alarm_compare

Mahalanobis threshold (baseline 99th percentile): 13.725

Mahalanobis alarm	False	True
PCA alarm
False	102	4
True	3	11

evaluation = []
for method in ["pca_alarm", "md_alarm"]:
    for cause, idx in anomaly_map.items():
        evaluation.append({
            "method": method,
            "cause": cause,
            "detected": int(df.loc[idx, method].sum()),
            "injected": len(idx),
        })
pd.DataFrame(evaluation)

	method	cause	detected	injected
0	pca_alarm	cooling_relation	2	4
1	pca_alarm	drive_relation	4	4
2	pca_alarm	pressure_relation	4	4
3	md_alarm	cooling_relation	4	4
4	md_alarm	drive_relation	4	4
5	md_alarm	pressure_relation	4	4

結果の読み取り

両手法が一致する点は優先度を上げ、不一致点は逸脱の性質を調べます。PCA再構成誤差は主要部分空間で説明できない関係崩れ、マハラノビス距離は中心からの総合的な遠さに反応します。検証用ラベルで検知件数を確認できますが、本番では故障記録や点検結果による継続的な評価が必要です。

No.069：多変量異常の原因候補を分析する

実務での意味

「異常です」だけでは点検につながりません。どのセンサーが通常構造から外れたかを示し、冷却系、駆動系、圧送系などの確認順序へ翻訳します。

分析・モデル化の考え方

PCAの標準化空間における残差

$$\mathbf{e}_i=\mathbf{z}_i-\hat{\mathbf{z}}_i$$

の二乗 $e_{ij}^2$ を変数別の寄与候補とします。これは因果寄与ではなく、再構成誤差をどの変数が構成したかを示す調査の手掛かりです。

Pythonで確認する

residual_sq = (X_all - X_reconstructed) ** 2
contribution = pd.DataFrame(residual_sq, columns=features, index=df.index)
df["priority_alarm"] = df["pca_alarm"] & df["md_alarm"]
priority_idx = df.index[(~train_mask) & df["priority_alarm"]]

root_cause_rows = []
for idx in priority_idx:
    top = contribution.loc[idx].sort_values(ascending=False).head(2)
    root_cause_rows.append({
        "timestamp": df.loc[idx, "timestamp"],
        "true_cause_for_validation": df.loc[idx, "true_cause"],
        "first_sensor_candidate": top.index[0],
        "second_sensor_candidate": top.index[1],
        "pca_error": df.loc[idx, "pca_error"],
        "mahalanobis_sq": df.loc[idx, "mahalanobis_sq"],
    })
root_cause_table = pd.DataFrame(root_cause_rows).sort_values("pca_error", ascending=False)
root_cause_table.head(12)

	timestamp	true_cause_for_validation	first_sensor_candidate	second_sensor_candidate	pca_error	mahalanobis_sq
4	2025-01-15 15:02:00	drive_relation	vibration_mm_s	coolant_flow_L_min	6.741	108.729
6	2025-01-15 15:04:00	drive_relation	vibration_mm_s	coolant_flow_L_min	6.655	122.452
3	2025-01-15 15:01:00	drive_relation	vibration_mm_s	coolant_flow_L_min	6.469	107.022
5	2025-01-15 15:03:00	drive_relation	vibration_mm_s	coolant_flow_L_min	3.355	60.361
9	2025-01-15 15:40:00	pressure_relation	pressure_MPa	coolant_flow_L_min	2.990	62.694
7	2025-01-15 15:38:00	pressure_relation	pressure_MPa	coolant_flow_L_min	2.671	72.714
10	2025-01-15 15:41:00	pressure_relation	pressure_MPa	temperature_C	2.274	66.094
8	2025-01-15 15:39:00	pressure_relation	pressure_MPa	coolant_flow_L_min	2.165	57.027
2	2025-01-15 14:28:00	cooling_relation	vibration_mm_s	coolant_flow_L_min	1.772	63.051
1	2025-01-15 14:26:00	cooling_relation	vibration_mm_s	motor_current_A	1.305	43.357
0	2025-01-15 14:14:00	normal	vibration_mm_s	pressure_MPa	0.836	15.219

top_idx = df.loc[(~train_mask) & df["priority_alarm"], "pca_error"].nlargest(8).index
plot_data = contribution.loc[top_idx]
fig, ax = plt.subplots(figsize=(9, 5))
bottom = np.zeros(len(plot_data))
for col in features:
    ax.bar(range(len(plot_data)), plot_data[col], bottom=bottom, label=col)
    bottom += plot_data[col].to_numpy()
ax.set_xticks(range(len(plot_data)), [df.loc[i, "timestamp"].strftime("%H:%M") for i in top_idx], rotation=45)
ax.set_title("Sensor-wise contribution candidates for priority alarms")
ax.set_xlabel("Alarm timestamp")
ax.set_ylabel("Squared standardized residual")
ax.grid(True, axis="y", alpha=0.3)
ax.legend(bbox_to_anchor=(1.02, 1), loc="upper left")
plt.tight_layout()
plt.show()

結果の読み取り

寄与の大きい2センサーを、設備図・制御ロジック・保全履歴と照合します。たとえば温度と冷却水流量が上位なら、熱交換器汚れ、バルブ開度、流量計ドリフトを候補にします。ただし、相関する変数間では寄与が分散するため、棒の最大値を故障原因と断定してはいけません。

No.070：高次元データにおける異常検知の注意点を整理する

実務での意味

センサー数を増やせば精度が自動的に上がるわけではありません。冗長センサー、欠損、時刻ずれ、モード混在は、距離や共分散推定を不安定にし、保守工数を増やします。

分析・モデル化の考え方

基準データ360件に対して、元の5変数へノイズ列を追加し、共分散行列の条件数とマハラノビス距離のばらつきを比較します。条件数が大きいほど逆行列計算は不安定になりやすく、データ件数に近い次元数では過学習リスクも高まります。

Pythonで確認する

dimension_check = []
for noise_dim in [0, 5, 20, 50, 100]:
    extra_train = rng.normal(size=(train_mask.sum(), noise_dim))
    extra_monitor = rng.normal(size=((~train_mask).sum(), noise_dim))
    high_train = np.column_stack([X_train, extra_train])
    high_monitor = np.column_stack([X_all[~train_mask], extra_monitor])
    high_cov = np.cov(high_train, rowvar=False)
    high_precision = np.linalg.pinv(high_cov)
    high_mu = high_train.mean(axis=0)
    centered_monitor = high_monitor - high_mu
    d2 = np.einsum("ij,jk,ik->i", centered_monitor, high_precision, centered_monitor)
    dimension_check.append({
        "variables": high_train.shape[1],
        "samples_per_variable": high_train.shape[0] / high_train.shape[1],
        "covariance_condition_number": np.linalg.cond(high_cov),
        "monitor_distance_cv": d2.std() / d2.mean(),
    })
dimension_check = pd.DataFrame(dimension_check)
dimension_check

	variables	samples_per_variable	covariance_condition_number	monitor_distance_cv
0	5	72.000	166.027	1.808
1	10	36.000	169.486	1.296
2	25	14.400	174.781	0.677
3	55	6.545	204.311	0.385
4	105	3.429	233.850	0.241

fig, ax1 = plt.subplots(figsize=(7, 4))
ax1.plot(dimension_check["variables"], dimension_check["samples_per_variable"], marker="o", color="#4C78A8")
ax1.set_title("Data support decreases as sensor dimensions grow")
ax1.set_xlabel("Number of variables")
ax1.set_ylabel("Baseline samples per variable")
ax1.set_yscale("log")
ax1.grid(True, which="both", alpha=0.3)
plt.tight_layout()
plt.show()

結果の読み取り

変数数が増えるほど1変数あたりの基準データが減り、共分散推定の信頼性が落ちます。実務では次を確認します。

• センサー時刻同期、欠損、校正履歴を整える
• 品種・工程・運転モードを分け、正常群を混ぜすぎない
• 目的に関係しない列や重複列を削減する
• PCA、縮小共分散、正則化などを候補にする
• 学習期間外の正常条件で誤報率を評価する
• 原因説明とモデル保守の担当を決める

この簡易実験の条件数は乱数に左右されるため、特定の数値を品質基準にするのではなく、次元増加に対する診断指標として扱います。

対象ノックを通して見える実務上の示唆

1. 単変量と多変量は併用する：安全限界はルール、多変量モデルは関係崩れの早期候補に使います。
2. 基準期間の品質がモデル品質を決める：正常と確認された期間、運転モード、保全直後などの情報が必要です。
3. 複数スコアで優先度を付ける：PCAとマハラノビス距離の一致を高優先度とし、不一致も捨てずに性質を調べます。
4. 説明は原因ではなく点検候補：変数別寄与を設備知識と結び付け、点検結果を記録します。
5. アラート件数をKPIにする：精度だけでなく、日次件数、確認時間、真因到達率、停止回避額を追います。

実務導入する場合に必要なこと

PoCから運用へ進む際は、次の設計が必要です。

領域	確認事項
データ	タグ定義、単位、周期、時刻同期、欠損、校正・交換履歴
基準期間	正常判定者、品種・負荷・季節の網羅、除外期間
評価	故障・品質記録との照合、見逃し・誤報コスト、期間外検証
アラート	しきい値、連続回数、抑制時間、優先度、通知先
業務	初動手順、点検チェックリスト、結果入力、責任分界
運用	分布変化監視、再学習条件、モデル版管理、監査ログ

まずは1設備・1故障モードに絞り、過去データで候補を再現できるか、1日のアラートを現場が処理できるかを確認する進め方が現実的です。

まとめ

No.061〜No.070では、相関と散布図による観察から始め、PCA再構成誤差とマハラノビス距離で多変量異常を数値化し、変数別残差で原因候補へつなげました。

重要なのは、異常スコアを故障診断と同一視しないことです。モデルは通常運転との差を示します。そこへ運転条件、設備構造、保全履歴、品質結果を重ねることで、初めて意思決定に使える点検優先度になります。

法人向けのご相談

数理工房では、製造業向けの異常検知PoC、センサーデータの品質診断、アラート設計、現場説明可能性を含む運用設計、企業研修をご支援しています。課題や保有データがまだ整理されていない段階でも、対象設備・意思決定・評価方法の切り分けからご相談いただけます。

📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。