問いと設定
異常な1件を説明することと、2期間の分布変化を説明することは異なるタスクです。DoWhy-GCMを使い、同じ5ノードのモデルに両方の問いを出します。
ambient→temp→pressure→vibrationの経路と、temp・pressure・vibrationからdefectへの経路を作ります。旧条件から新条件への変更は、ambientの平均を1上げることと、pressure←tempの係数を0.8から1.2へ変えることです。実在の設備やデータではありません。
各ノードに線形の加法ノイズ機構を置きます。構造は既知で、外生ノイズは独立です。新旧データは同じ乱数から生成して比較の揺れを抑えていますが、実際の時系列監視データではないことに注意してください。
実行できるコード
Python 3.12.14で実行しました。解析はすべて合成データを使う研究・学習用の例です。
専用の仮想環境で、次の版をインストールして実行します。数値は2026-10-02の実行結果です。
python -m pip install numpy==2.4.6 pandas==3.0.6 scipy==1.15.3 scikit-learn==1.6.1 statsmodels==0.15.0 dowhy==0.14
次のコードは上から一括実行できます。seedを固定していますが、依存ライブラリや実装の変更で細かな数値は変わり得ます。
import numpy as np
import pandas as pd
import networkx as nx
from dowhy import gcm
from dowhy.gcm.shapley import ShapleyConfig
def generate(shift=False):
rng = np.random.default_rng(20261002)
n = 3000
ambient, et, ep, ev, ed = rng.normal(size=(5,n))
ambient += float(shift)
temp = .9*ambient+.4*et
pressure = (.8+.4*shift)*temp+.4*ep
vibration = .7*pressure+.4*ev
defect = .6*temp+.9*pressure+1.1*vibration+.4*ed
return pd.DataFrame(dict(ambient=ambient,temp=temp,pressure=pressure,vibration=vibration,defect=defect))
edges = [('ambient','temp'),('temp','pressure'),('pressure','vibration'),
('temp','defect'),('pressure','defect'),('vibration','defect')]
def fit(data):
model = gcm.InvertibleStructuralCausalModel(nx.DiGraph(edges))
model.set_causal_mechanism('ambient',gcm.EmpiricalDistribution())
for node in ['temp','pressure','vibration','defect']:
model.set_causal_mechanism(node,gcm.AdditiveNoiseModel(gcm.ml.create_linear_regressor()))
gcm.fit(model,data)
return model
old,new = generate(),generate(True)
old_model,new_model = fit(old),fit(new)
np.random.seed(20261002)
anomaly = pd.DataFrame(dict(ambient=[.2],temp=[.25],pressure=[2.7],vibration=[2.2],defect=[6.]))
scores = gcm.attribute_anomalies(old_model,'defect',anomaly,
anomaly_scorer=gcm.MedianCDFQuantileScorer(),num_distribution_samples=1000,
shapley_config=ShapleyConfig(n_jobs=1))
print('Anomaly score contributions:',{k:round(float(np.asarray(v).reshape(-1)[0]),4) for k,v in scores.items()})
np.random.seed(20261002)
change = gcm.distribution_change_of_graphs(old_model,new_model,'defect',num_samples=2000,
difference_estimation_func=lambda x,y:float(np.mean(y)-np.mean(x)),
shapley_config=ShapleyConfig(n_jobs=1))
print('Observed mean difference:',round(new.defect.mean()-old.defect.mean(),4))
print('Signed mean-change contributions:',{k:round(float(v),4) for k,v in change.items()})
print('Sum:',round(sum(change.values()),4))
実行結果
| ノード | 異常スコアへの寄与 | 平均分布変化への寄与 |
|---|---|---|
| ambient | −1.4485 | 2.0346 |
| temp | −0.2365 | 0.0020 |
| pressure | 5.0270 | 0.2933 |
| vibration | 0.3358 | −0.0314 |
| defect | 1.3587 | −0.0541 |
異常列は、指定した1行のdefectに対するMedianCDFQuantileScorerの寄与です。分布変化列は旧→新のoutcome平均の符号付き変化を分解した値で、単位も尺度も異なります。列同士の大きさは比較しません。
新旧標本の実際のdefect平均差は2.3388、モデルからMonte Carloで分解した寄与の合計は2.2444でした。差は約0.0944あり、標本の平均差を厳密に再現した分解とは書けません。
分解の対象を先に決める
異常要因解析は、観測行から復元した外生ノイズが異常スコアにどう寄与するかを扱います。pressureの寄与が最大でも、実設備のpressureに障害があると立証した結果ではありません。指定したDAG、学習機構、スコアの下での説明です。
分布変化の分解は、旧機構と新機構を入れ替えたモデルを比較し、Shapley値として寄与を配分します。今回は平均差を符号付きで指定しました。絶対平均差を使うと別の量を分解するため、結果を同じ意味で読めません。
負の寄与と残る誤差
Shapley寄与は負にもなります。指定したスコアや統計量を下げる方向への配分を表し、負の故障確率という意味ではありません。平均分布変化の小さい負値には、学習誤差・Monte Carlo誤差も含まれます。
今回の寄与は点推定です。反復seed、生成標本数、学習標本、モデル仕様を変えた安定性や不確実性を評価していないので、細かな順位を確定しません。既知の変更先ambient・pressureが主な寄与として現れることを確認した教材です。
現場へ適用する場合には、同じ対象・測定方式の新旧データか、選択や欠測が変わっていないか、機構モデルの残差が妥当かを確認します。異常診断と分布ドリフトを混ぜず、説明したい対象と統計量を記録する必要があります。