体温を測って「37.2℃」という数値を見せられても、それだけでは平熱か発熱か判断できません。平熱域は個人差があり、朝と夕方でも変わるからです。同じように、あるセンサ値が「異常かどうか」は、それ単体では決まらないことがよくあります。その値が観測された「状況」「文脈」を知って初めて、異常かどうかが決まるのです。
この発想を初めて厳密な統計モデルとして定式化したのが、本記事で扱う Xiuyao Song, Mingxi Wu, Christopher Jermaine, Sanjay Ranka, “Conditional Anomaly Detection,” IEEE Transactions on Knowledge and Data Engineering, vol. 19, no. 5, pp. 631–645, 2007 です。原論文を精読し、数式・実験プロトコル・実験結果を可能な限り原文に忠実に、しかし直感的にたどり直します。
なぜこれを学ぶのか
論文の Abstract は、実務家にとって示唆に富む問題提起から始まります。異常検知ソフトウェアがデータ分析ツールとして使われる場合、最も検知しづらい異常を見つけることが最重要課題とは限らない、というのです。むしろ重要なのは、システムが「異常」として報告してくる候補が、実際にユーザーの関心を引く、興味深いものであることです。あまりに凡庸な点まで異常候補として返してしまうと、そのソフトウェアはすぐに使われなくなってしまいます。
この問題の根っこにあるのが、属性を区別しないという暗黙の仮定です。統計的な異常検知手法の多くは「全データが同一の生成過程から作られた」という帰無仮説を検定する形を取りますが、その際「一部の属性はそもそも異常の直接的な証拠として扱うべきではない」という事前知識を活かせていません。応用先は2つ考えられます。
- 公衆衛生サーベイランス：疾病のアウトブレイクをできるだけ早く検知したい。気温のような環境変数と、発熱者数のような指標変数を区別しないと、単に「暑い日/寒い日」というだけで大量の誤警報が出てしまいます。
- 一般的な異常検知システムの運用：人間の担当者が確認できる「異常候補の予算」には限りがあります。環境的に珍しいだけの点まで候補に含めると、その予算が無駄遣いされ、本当に注意すべき異常が埋もれてしまいます。
この記事の内容
- 環境属性・指標属性の区別と、原論文 Fig.1 の具体例（気温と発熱者数）
- 生成モデル $f_{CAD}$ の3ステップの物語と、その数式
- 異常判定のカットオフ（予算 $\epsilon$）の仕組み
- EMアルゴリズムによる学習 — 責任度 $b_{kij}$ の導出と更新式
- 3つの学習アルゴリズム（Direct-CAD／GMM-CAD-Full／GMM-CAD-Split）の正確な違い
- 論文の摂動評価プロトコルと、13個の実データセットでの実験結果
- Python実装：対応づけの不動点反復EMをscikit-learnで実装し、原論文 Fig.1 のシナリオを数値で再現検証
前提知識
1. 環境属性と指標属性 — 原論文 Fig.1 の具体例
論文はまず、感染症サーベイランスの具体例を使って動機を説明します。2つの変数を毎日観測しているとしましょう。
- max_daily_temp：その日の最高気温
- num_fever：その日、発熱を訴えて救急外来に来た人数
出典: Song, Wu, Jermaine, Ranka, “Conditional Anomaly Detection,” IEEE TKDE, 2007, Fig. 1.
原論文はこの散布図から2つの点を対比させます。点Aは最高気温が「very cold」の極端な側にあり、発熱者数も「many」（最大級）です。単純な外れ値検知（例えば全体に対して低密度な点を探す方法）であれば、点Aは間違いなく最も明白な外れ値として検出されるでしょう。点Bは最高気温が中庸な範囲にありながら、発熱者数が周囲の点よりわずかに高い位置にあります。
ここでの論文の主張が核心です。「もし max_daily_temp が異常の直接的な証拠ではないという追加情報を活用できるなら、点Aを無視するのはむしろ安全である」。なぜか。ある June の暑い日に気温が70°Fしかなかった（平均は90°F近い）としても、それはただの「異常に涼しい日」であって、アウトブレイクの証拠にはなりません。逆に、寒い日に発熱者数が多いのは、寒さで免疫力が落ちる、あるいは単に季節性のかぜが増えるという、よく知られた自然な相関に過ぎない可能性が高いのです。一方、点Bのように「気温からは説明できない発熱者数の増加」こそ、真に注意すべきシグナルです。
再現用に、この状況を自分で作った連続データで可視化してみます。
破線は気温から期待される発熱者数の傾向（寒くても暑くても発熱者は少し増えるU字カーブ）です。点Aは気温という説明変数の値こそ極端ですが、傾向線からのズレは小さく、条件付きに見れば正常です。点Bは傾向線から明確に浮いており、これが本当に調べるべき逸脱です。「指標属性の値そのものの珍しさ」ではなく「環境属性から期待される値からのズレ」を測る——これが条件付き異常検知のすべての出発点です。
論文はこの区別に基づき、データ属性を2種類に分けます。
- 環境属性（environmental attributes）$x$：運用モード・操作・気象条件など。それ自体を異常の直接的な証拠として扱わないという事前知識をユーザーが与える量。
- 指標属性（indicator attributes）$y$：実際に異常かどうかを判定したい量。
重要なのは、環境属性を「無視する」わけではないという点です。環境属性は「指標属性の期待値にどう影響するか」を通じてのみモデルに効いてきます。次章で、この考えを厳密な生成モデルとして組み立てます。
2. 生成モデル $f_{CAD}$ — 3ステップのサイコロ
論文は $f_{CAD}(y\mid\Theta, x)$ という条件付き密度関数を、次の3ステップの生成過程として定義します。
出典: Song, Wu, Jermaine, Ranka, “Conditional Anomaly Detection,” IEEE TKDE, 2007, Fig. 2.
図の左側が環境属性の空間、右側が指標属性の空間です。左には環境属性の混合ガウスモデル $U$ の3成分 $U_1, U_2, U_3$ が等高線で描かれ、右には指標属性の混合ガウスモデル $V$ の3成分 $V_1, V_2, V_3$ が描かれています。ある環境属性の値 $x$（黒丸）が観測されたとき、対応する指標属性 $y$（黒丸）が生成される過程は次の3段階です。
- Step 1：$x$ を生成したのがどのガウス成分 $U_i$ かを、事後確率 $p(x \in U_i)$ に基づいて（仮想的に）決める。図の例では $U_2$ が選ばれています。
- Step 2：$U_2$ から $V$ 側のどの成分に「マッピングされるか」を、確率 $p(V_j \mid U_2)$ に従ったサイコロで決める。図の例では $p(V_1\mid U_2)=.2,\ p(V_2\mid U_2)=.2,\ p(V_3\mid U_2)=.6$ であり、$V_3$ が選ばれています。
- Step 3：選ばれた $V_3$ から $y$ を1点サンプリングする。
この物語で重要なのは、この生成過程は $x$ を実際に生成するわけではないという点です。$x$ は既に観測された値として与えられ、あくまで「$x$ を生成したのがどの $U_i$ か」という推論のためだけに $U$ が使われます。$U$ の役割は「$x$ の分布を説明すること」ではなく「$x$ を手がかりにして、どの $V$ 成分が期待される指標分布かを決めるための橋渡し」です。
この3ステップを数式にすると、条件付き密度は次のようになります。
$$ \begin{equation} f_{CAD}(y\mid\Theta, x) = \sum_{i=1}^{n_U} p(x \in U_i) \sum_{j=1}^{n_V} f_G(y\mid V_j)\, p(V_j\mid U_i) \end{equation} $$
ここで、
- $f_G(y\mid V_j)$ は指標属性側の $j$ 番目のガウス成分 $V_j = \langle \mu_{V_j}, \Sigma_{V_j}\rangle$ の確率密度関数
- $p(V_j\mid U_i)$ は「$x$ が $U_i$ から生成されたと分かっているとき、指標属性が $V_j$ から生成される確率」——モデルパラメータの一部として直接学習される対応づけ
- $p(x \in U_i)$ は「$x$ が環境側の $i$ 番目のガウス成分から生成された事後確率」で、ベイズの定理から計算できます
$$ \begin{equation} p(x \in U_i) = \frac{f_G(x\mid U_i)\, p(U_i)}{\sum_{k=1}^{n_U} f_G(x\mid U_k)\, p(U_k)} \end{equation} $$
式を読み解くと、Σが二重になっている理由がはっきりします。外側の $\sum_i$ は「$x$ がどの環境成分から来た可能性があるか」を重み $p(x\in U_i)$ で加重平均し、内側の $\sum_j$ は「その環境成分から見て、指標属性がどの $V_j$ に対応しやすいか」を $p(V_j\mid U_i)$ で加重してガウス密度を足し合わせる、という環境 → 対応する指標 → その尤度という2段階の橋渡しです。パラメータの全体集合は $\Theta = \{U, V, p(V_j\mid U_i)\}$ の3つ組になります。
この生成モデルの美しさは、環境属性の「珍しさ」自体は $f_{CAD}$ に直接効かないという点にあります。$p(x\in U_i)$ はあくまで「複数の $U_i$ 候補の中でどれが相対的にもっともらしいか」という相対的な重みであり、絶対的な密度 $f_G(x\mid U_i)$ の大小そのものは最終的な $f_{CAD}$ の値には現れません。$x$ がどれほど珍しくても、対応する $V_j$ から見て $y$ がもっともらしければ、$f_{CAD}$ は高い値を取ります。これがまさに、点Aを安全に無視できる仕組みの数学的な正体です。
3. 異常の判定 — 尤度による順位づけとカットオフ
$f_{CAD}$ が定まれば、異常検知の考え方は単純です。パラメータ $\Theta$ が最適化されていれば、$f_{CAD}$ は過去から未来までのあらゆるデータ点に対して「驚きの少なさ」の意味のある順位づけを与えます。小さい値を取る点ほど、データ空間の低密度領域にあり、異常として疑うべきです。
ただし、実務では「何を下回ったら異常と呼ぶか」というカットオフ値を決める必要があります。論文はここで、冒頭の「異常候補の予算」という発想に立ち返ります。ユーザーはまず、0から1の間の割合 $\epsilon$ を選びます。$\epsilon$ が小さいほど、検知は選択的（厳しい）になります。次に、訓練データの $f_{CAD}$ 値を昇順にソートし、$c = \epsilon n$ 番目（$n$ は訓練データ数）の点をカットオフ地点とします。新しい観測点 $(x_{new}, y_{new})$ が来たとき、
$$ \begin{equation} f_{CAD}(y_{new}\mid\Theta, x_{new}) < f_{CAD}(y_c\mid\Theta, x_c) \end{equation} $$
が成り立てば異常フラグを立てます。これは「訓練データの中で $\epsilon$ 割よりもさらに驚くべき値が観測されたら知らせる」という、実務家にとって直感的な運用そのものです。$\epsilon$ を大きくすれば見逃しは減りますが誤報が増え、小さくすればその逆になります。
4. パラメータ学習 — EMアルゴリズムと責任度 $b_{kij}$
生成モデル $f_{CAD}$ を定義しただけでは不十分で、これを既存の（正常な）データにフィットさせる手続きが必要です。論文はここで 最尤推定（MLE） の枠組みを使います。目的関数は対数尤度
$$ \begin{equation} \Lambda = \sum_{k=1}^{n} \log f_{CAD}(y_k\mid\Theta, x_k) \end{equation} $$
で、これを $\Theta$ について最大化します。しかし $f_{CAD}$ の中には「$x_k$ を生成した $U_i$」「$y_k$ を生成した $V_j$」という観測されない（隠れた）変数が埋め込まれているため、直接の最大化は困難です。ここで EMアルゴリズム（Dempster, Laird, Rubin, 1977）の出番です。
EMは「隠れ変数の真の値」を知る代わりに、現在のパラメータの下での隠れ変数の事後分布に関して目的関数の期待値を取り（E-step）、その期待値を最大化するようにパラメータを更新する（M-step）ことを繰り返します。この手続きは大域最適を保証しませんが、各反復で尤度が単調に改善することは保証されます。
条件付き異常検知の場合、隠れ変数は「$x_k$ を生成した $U$ の成分 $i$」と「$y_k$ を生成した $V$ の成分 $j$」の組 $(i,j)$ です。この同時事後確率を論文は $b_{kij}$ と表記します。
$$ \begin{equation} b_{kij} = \frac{f_G(x_k\mid U_i)\, p(U_i)\, f_G(y_k\mid V_j)\, p(V_j\mid U_i, \Theta)} {\displaystyle\sum_{t=1}^{n_U}\sum_{h=1}^{n_V} f_G(x_k\mid U_t)\, p(U_t)\, f_G(y_k\mid V_h)\, p(V_h\mid U_t,\Theta)} \end{equation} $$
読み方は「データ点 $k$ について、環境成分 $i$ と指標成分 $j$ の組が同時に真である事後確率」で、分子はその組の（環境側×指標側の）尤度、分母はあらゆる $(i,j)$ の組についての正規化項です。ここで $p(V_j\mid U_i,\Theta)$ が現在の推定値であることに注意してください——これはEMの各反復で更新される対応づけそのものであり、$b_{kij}$ の計算はこの現在値に依存します。
M-stepでは、この $b_{kij}$ を重みとして、通常のGMMのEMと同じ「加重平均・加重共分散」の形の更新式が導けます。ラグランジュ未定乗数法を使い、Q関数を4つの独立な項に分解して1つずつ最大化するという手順です（詳細な導出は原論文 Appendix B に譲ります）。結果として得られる更新式は次の通りです。
$$ \begin{align} \overline{p(U_i)} &= \frac{\sum_k \sum_j b_{kij}}{\sum_k \sum_h \sum_j b_{khj}} \ , \qquad \bar\mu_{U_i} = \frac{\sum_k \sum_j b_{kij}\, x_k}{\sum_k \sum_j b_{kij}} \ , \qquad \bar\Sigma_{U_i} = \frac{\sum_k \sum_j b_{kij}(x_k-\bar\mu_{U_i})(x_k-\bar\mu_{U_i})^T}{\sum_k \sum_j b_{kij}} \\[4pt] \bar\mu_{V_j} &= \frac{\sum_k \sum_i b_{kij}\, y_k}{\sum_k \sum_i b_{kij}} \ , \qquad \bar\Sigma_{V_j} = \frac{\sum_k \sum_i b_{kij}(y_k-\bar\mu_{V_j})(y_k-\bar\mu_{V_j})^T}{\sum_k \sum_i b_{kij}} \ , \qquad \overline{p(V_j\mid U_i)} = \frac{\sum_k b_{kij}}{\sum_k \sum_h b_{kih}} \end{align} $$
式変形の骨子だけ捉えると、$U_i$ のパラメータは「$j$ について周辺化した責任度 $\sum_j b_{kij}$」で重みづけた $x_k$ の加重平均・加重共分散であり、$V_j$ のパラメータは逆に「$i$ について周辺化した責任度 $\sum_i b_{kij}$」で重みづけた $y_k$ の加重平均・加重共分散です。そして対応づけ $p(V_j\mid U_i)$ は、$U_i$ に紐づく責任度の総量のうち $V_j$ に流れた割合として計算されます。この最後の式は、通常のGMMのEMには存在しない、条件付き異常検知に固有の更新式です。
5. 3つの学習アルゴリズム — 何を「同時に」学ぶかが違う
ここまでの定式化を素直に実装すると、$U$・$V$・対応づけ $p(V_j\mid U_i)$ の3つを1本のEMで同時に最適化することになります。これが論文の第一の手法 Direct-CAD です。しかし、$U$ と $V$ という2組のガウス混合と、それらを結ぶ写像を同時に最適化する問題は非常に複雑で、EMは局所最適解に収束しやすくなります。
そこで論文は、計算をより単純な部分問題に分割する2つの代替アルゴリズムを提案します。3つの違いを図にすると次のようになります。
- Direct-CAD：前述の通り、環境データと指標データから $U$・$V$・対応づけを1本のEMで同時最適化します。最も理論的に忠実ですが、計算コストが高く局所解に陥りやすいという欠点があります。
- GMM-CAD-Full：まず環境属性と指標属性を連結したベクトル $(x,y)$ に対して、$(d_U + d_V)$ 次元の単一の混合ガウスモデルを学習します（成分数は $n_U = n_V$ と仮定）。学習された各ガウス成分の平均・共分散を、環境属性の次元だけ取り出せば $U$ に、指標属性の次元だけ取り出せば $V$ になります（単なる射影で、追加の計算は不要です）。$U$ と $V$ が定まった後、対応づけ $p(V_j\mid U_i)$ だけを学習する簡易版のEM（4節の $b_{kij}$ と最後の更新式のみを反復）を走らせます。この方法は「$U$ と $V$ の各成分が共分散構造を通じて自然に対応づけられている」という直感に基づいています——同じ結合ガウス成分から射影されたのだから、対応するのは当然だ、というわけです。
- GMM-CAD-Split：$U$ と $V$ を完全に独立な2本のEMで別々に学習します（環境データだけでGMMを学習し、指標データだけで別のGMMを学習する）。そのうえで、GMM-CAD-Fullと同じ「対応づけだけを学習する簡易EM」を走らせます。3つの中で最も軽量な計算量で済みます。
論文の実験結果（次章で詳しく見ます）では、GMM-CAD-Full が総合的に最も高い精度を示しました。直感的には、結合空間で学習することで環境と指標の共分散構造がある程度自然に反映されるためだと考えられますが、論文自身も「どのデータセットでも一様に他手法より優れているわけではない」と正直に注記しています。
計算量の面でも3手法は異なります。Direct-CADは1反復あたり $O(n\, n_U n_V (d_U^2+d_V^2))$、GMM-CAD-Fullは結合GMMの学習に $O(n\, n_U (d_U+d_V)^2)$、対応づけの学習に $O(n\, n_U^2)$、GMM-CAD-Splitは環境側 $O(n\, n_U d_U^2)$・指標側 $O(n\, n_V d_V^2)$ を別々に、という具合です。本記事のPython実装では、最も理解しやすく実装が軽いGMM-CAD-Splitを採用します。
6. 実験プロトコル — 「摂動」で条件付き異常を人工的に作る
条件付き異常検知の評価には工夫が要ります。実データには「これは条件付き異常である」というラベルがほとんど存在しないためです。論文が考案した評価プロトコルは、この論文のもう一つの重要な貢献です。
各データセットについて、次の手順を10回繰り返します。
- データの80%を訓練データ、20%をテストデータ $testData$ とする。
- $testData$ のうち20%を、環境属性の値だけを見た標準的なGMMベースの外れ値検定でランキングし、上位（最も低密度な）ものを $outliers$ と呼ぶ。ここで指標属性は一切無視される点が重要です。
- $testData$ を $perturbed$ と $nonPerturbed$ の2つの等サイズな集合に分割する。ただし $outliers \subset nonPerturbed$ となるよう強制する。
- $perturbed$ の各点について、指標属性の値を入れ替える。
- $testData = perturbed \cup nonPerturbed$ として異常検知を実行する。
ここでの入れ替え（摂動）の手順が巧妙です。単純にノイズを加えるのではなく、各点 $z=(x,y) \in D$ について、$k=\min(50, |D|/4)$ 個の点を無作為にサンプリングし、その中で指標属性がユークリッド距離で最も遠い点 $z’=(x’, y’)$ を選び、新しい点 $(x, y’)$ を作ります。つまり環境属性は本物のまま、指標属性だけを「その環境ではまず起こらないであろう」別の値にすり替えるわけです。これはまさに条件付き異常の定義そのもの——環境属性 $x$ とその環境下では期待されない指標属性 $y’$ の組み合わせ——を人工的に作る操作になっています。
評価の理想形は明確です。
- $perturbed$ の点は高い割合で異常と判定されてほしい（真の条件付き異常だから）
- $nonPerturbed$ の点は低い割合で異常と判定されてほしい（元の分布からのサンプルのままだから）
- $outliers$（環境属性だけが珍しい点。$nonPerturbed$ の部分集合）も同じくらい低い割合で異常と判定されてほしい——これが最大のポイントで、環境属性の珍しさだけに釣られて誤検知を出す手法をあぶり出す設計になっています
7. 実データセットでの実験結果
論文はこのプロトコルを、次の13個の実データセット（環境属性数×指標属性数）に適用しています。
Synthetic (50×50, CADの生成モデルから合成)、Algae (11×6, UCI KDD Archive の藻類データ)、Streamflow (205×100, 米国の河川流量と気象データ)、ElNino (4×5, 太平洋の海洋観測ブイデータ)、Physics (669×70, KDD Cup 2003 の物理論文キーワード頻度)、Bodyfat (13×2, CMU StatLibの体脂肪率)、Houses (8×1, カリフォルニアの住宅価格)、Boston (15×1, ボストンの住宅価格)、FCAT-math / FCAT-reading (14×12, 14×11, フロリダの学力テストと学校特性)、FLFarms / CAFarms (114×52, 115×51, フロリダ/カリフォルニアの農場市場価値)、CApeaks (2×1, カリフォルニアの山岳標高と位置)。
比較対象は、単純なGMM、$k$-NN（$k=5$、最近傍距離ベースの外れ値スコア）、LOF（Local Outlier Factor）、そして本記事で見た3つのCAD手法（Direct-CAD／GMM-CAD-Full／GMM-CAD-Split）の計6手法です。
Table 1（$perturbed$ を正しく異常と判定できたか、13データセット平均の recall/precision）では、GMM-CAD-Full が平均 0.793 で最良でした（GMM 0.720、CAD-Split 0.737、Direct CAD 0.730、5th-NN 0.721、LOF 0.721）。ヘッド・トゥ・ヘッドの比較でも、GMM-CAD-Full は他の全手法に対して13データセット中10〜11勝という結果を残しています。
Table 2（$outliers$ を正しく非異常と判定できたか）では差がさらに開きます。GMM-CAD-Full が平均 0.749 で最良、対して単純GMMと5th-NNは平均 0.500（ランダムな判定と変わらない！）にとどまりました。これはまさに、「環境属性を条件づけずに扱う手法は、環境が珍しいだけの点を無差別に異常と誤検知してしまう」という論文の問題提起を裏付ける結果です。特に1,000次元近い Physics データセットでは、LOFと5th-NNの precision/recall がランダムなラベリングより悪くなる（50%未満）という顕著な結果も報告されています。
論文はさらに Appendix A で、ラベル付きの KDD Cup 1999 データセットに対する10分割交差検証も行っています。
出典: Song, Wu, Jermaine, Ranka, “Conditional Anomaly Detection,” IEEE TKDE, 2007, Fig. 3.
このデータセットに対しては LOF 以外の全手法がほぼ同等の高い precision（0.85〜1.0近辺）を示しており、論文自身も「このタスクは比較的簡単で、追加の結論を引き出すのは難しい」と述べています。実世界の侵入検知データでは属性間の相互作用がそれほど複雑ではなかった、というのが妥当な解釈でしょう。
総合すると、論文の結論は「どの手法も万能ではないが、任意のデータセットに対して最も期待できる選択肢を1つ挙げるなら GMM-CAD-Full である」というものです。
8. Python実装 — GMM-CAD-Splitと対応づけの不動点反復EM
理論を実装で確かめます。環境属性 $x$（1次元）・指標属性 $y$（2次元）を持つ合成データを作り、真の対応が「環境クラスタ $i$ → 指標クラスタ $i$」であるようにします。学習にはGMM-CAD-Split（$U$、$V$を独立に学習し、対応づけだけを別途学習する最も軽量な方式）を使います。
import numpy as np
from sklearn.mixture import GaussianMixture
from scipy.stats import multivariate_normal as mvn
rng = np.random.default_rng(0)
KU, KV = 3, 3
x_means = np.array([-3.0, 0.0, 3.0]); x_sd = 0.5
y_means = np.array([[-2.5, -2.0], [0.5, 2.6], [2.8, -1.6]])
y_covs = [(lambda A: A@A.T*0.3 + 0.15*np.eye(2))(rng.normal(0,1,(2,2))) for _ in range(KV)]
N = 2400
zs = rng.integers(KU, size=N)
x_all = rng.normal(x_means[zs], x_sd).reshape(-1, 1)
y_all = np.array([rng.multivariate_normal(y_means[z], y_covs[z]) for z in zs])
n_tr = 1700
x_tr, y_tr = x_all[:n_tr], y_all[:n_tr]
x_te_n, y_te_n, z_te_n = x_all[n_tr:], y_all[n_tr:], zs[n_tr:]
テスト側の「摂動」データは、論文のプロトコルを単純化し、環境はそのまま・指標だけ別クラスタ寄りにずらして作ります。
na = len(x_te_n)
y_an = []
for i in range(na):
zi = z_te_n[i]; w = (zi + 1) % KU
tgt = y_means[zi] + 0.85 * (y_means[w] - y_means[zi])
y_an.append(rng.multivariate_normal(tgt, y_covs[w]))
y_an = np.array(y_an); x_an = x_te_n.copy()
x_te = np.vstack([x_te_n, x_an]); y_te = np.vstack([y_te_n, y_an])
label = np.concatenate([np.zeros(na), np.ones(na)]) # 1=条件付き異常
$U$・$V$ を独立に学習した後、対応づけ $p(V_j\mid U_i)$ だけを学習する簡易EM（4節の $b_{kij}$ と更新式を、$U$・$V$ 自体は固定したまま繰り返す不動点反復）を実装します。ここが旧版の実装で単純化しすぎていた部分で、論文の式(34)が本来求めているのは1回の共起計算ではなく、この反復による収束解です。
gU = GaussianMixture(KU, covariance_type="full", random_state=0).fit(x_tr)
gV = GaussianMixture(KV, covariance_type="full", random_state=0).fit(y_tr)
pU_tr = gU.predict_proba(x_tr) # p(x_k∈U_i)、固定
dens_tr = np.stack([mvn.pdf(y_tr, gV.means_[j], gV.covariances_[j])
for j in range(KV)], axis=1) # f_G(y_k|V_j)、固定
M = np.full((KU, KV), 1.0 / KV) # p(V_j|U_i) の初期値(一様)
for it in range(50):
joint = pU_tr[:, :, None] * dens_tr[:, None, :] * M[None, :, :] # b_kij(未正規化)
b = joint / joint.sum(axis=(1, 2), keepdims=True)
M_new = b.sum(axis=0) # Σ_k b_kij
M_new = M_new / M_new.sum(axis=1, keepdims=True) # 式(34): 行(iごと)に正規化
if np.abs(M_new - M).max() < 1e-8:
M = M_new; break
M = M_new
print("収束した対応づけ p(V_j|U_i):\n", np.round(M, 3), " /", it+1, "反復")
実行結果：
収束した対応づけ p(V_j|U_i):
[[1. 0. 0.]
[0. 0. 1.]
[0. 1. 0.]] / 5 反復
学習された対応づけ行列はきれいな置換行列（対角の並べ替え）になっています。真の対応「環境 $i$ → 指標 $i$」を、教師なしのまま正しく復元できていることが分かります。GMMの成分番号の付与順は学習ごとに変わるため、対角そのものではなく置換行列になっていますが、構造としては完全に正しく学習されています。
不動点反復の収束の様子も見ておきましょう。
反復1回目から2回目で対数尤度が大きく改善し、その後は5回でぴたりと収束しています。$U$・$V$ をあらかじめ固定した簡易版のEMなので、フルの同時最適化より収束が速いことが分かります。これがGMM-CAD-SplitやGMM-CAD-Fullが計算量の面で有利な理由です。
続いて $f_{CAD}$ を実装し、無条件モデル（指標属性のみを見たGMM）と比較します。
def f_cad(xq, yq):
pUq = gU.predict_proba(xq)
densq = np.stack([np.atleast_1d(mvn.pdf(yq, gV.means_[j], gV.covariances_[j]))
for j in range(KV)], axis=1)
return ((pUq @ M) * densq).sum(axis=1) # Σ_i p(x∈U_i) Σ_j p(V_j|U_i) f_G(y|V_j)
s_cond = -np.log(f_cad(x_te, y_te) + 1e-300)
s_uncond = -gV.score_samples(y_te)
def auc(score, y):
order = np.argsort(score); ranks = np.empty(len(score)); ranks[order] = np.arange(1, len(score)+1)
npos = y.sum(); nneg = len(y) - npos
return (ranks[y == 1].sum() - npos*(npos+1)/2) / (npos*nneg)
print(f"無条件(指標yのみ) AUC={auc(s_uncond, label):.3f}")
print(f"条件付き(GMM-CAD-Split) AUC={auc(s_cond, label):.3f}")
無条件(指標yのみ) AUC=0.726
条件付き(GMM-CAD-Split) AUC=1.000
指標属性だけを見た無条件モデルでは、摂動データの指標値がたまたま別の正常クラスタと重なってしまうケースを見逃し、AUCは0.726にとどまります。一方、環境属性で条件づけたモデルは「この環境クラスタからは、この指標値は出るはずがない」を正確に捉え、ほぼ完全な分離（AUC=1.000）を達成しています。
密度そのものを可視化すると、この違いがより直感的に分かります。
無条件密度 $p(y)$（左）では3つの指標クラスタすべてが等しく「正常」として現れるため、別の環境クラスタに本来対応する指標値でもそれなりの密度を得てしまいます。条件付き密度 $p(y\mid x=\text{環境クラスタ1})$（右）では、対応するクラスタだけに密度の山が立ち、それ以外の領域はほぼゼロになります。環境で条件づけることで「正常領域」が大幅に絞り込まれるのです。
9. 原論文 Fig.1 のシナリオを数値で再現する
最後に、1節で見た「点A（環境が極端でも指標は期待通り）」と「点B（環境は平均的なのに指標が不整合）」の状況を、学習済みモデルで実際にスコアリングして検証します。ここでは $f_{CAD}$ に加えて、環境属性と指標属性を区別せず同列に扱う素朴な結合GMMも比較対象に加えます。これは論文が問題提起した「すべての属性を平等に扱う既存手法」に相当するベースラインです。
from sklearn.mixture import GaussianMixture as GM
g_joint = GM(KU, covariance_type="full", random_state=0).fit(np.hstack([x_tr, y_tr]))
# 点A: 環境が訓練分布の裾(極端に寒い)だが、指標はそのクラスタの期待通り
xA = np.array([[x_means[0] - 3.2 * x_sd]])
yA = (y_means[0] + rng.multivariate_normal([0, 0], y_covs[0]*0.15)).reshape(1, -1)
# 点B: 環境は平均的(クラスタ1)だが、指標は別クラスタ(クラスタ2)寄り
xB = np.array([[x_means[1]]])
yB = (y_means[2] + rng.multivariate_normal([0, 0], y_covs[2]*0.15)).reshape(1, -1)
all_fcad = f_cad(x_all, y_all)
all_joint = np.exp(g_joint.score_samples(np.hstack([x_all, y_all])))
all_v = np.exp(gV.score_samples(y_all))
def pct(val, arr): # 正常データ中でのパーセンタイル(低いほど異常寄り)
return (arr < val).mean() * 100
for name, xq, yq in [("点A", xA, yA), ("点B", xB, yB)]:
fc = f_cad(xq, yq)[0]
fj = np.exp(g_joint.score_samples(np.hstack([xq, yq]))[0])
fv = np.exp(gV.score_samples(yq)[0])
print(name, "f_CAD%ile=", round(pct(fc, all_fcad), 1),
" joint%ile=", round(pct(fj, all_joint), 1),
" f_V%ile=", round(pct(fv, all_v), 1))
点A f_CAD%ile= 99.0 joint%ile= 1.7 f_V%ile= 99.0
点B f_CAD%ile= 0.0 joint%ile= 0.0 f_V%ile= 83.0
数値は原論文の主張を見事に裏付けています。
- 点A（環境は極端だが指標は期待通り）：条件付き $f_{CAD}$ は99.0パーセンタイル（極めて正常）と正しく判定します。ところが環境と指標を同列に扱う結合GMMは1.7パーセンタイルという最下位クラスの値を出し、誤って異常フラグを立ててしまいます。環境属性の値そのものが訓練データの裾にあるため、結合空間ではこの点全体が低密度領域に入ってしまうのです。
- 点B（環境は平均的だが指標が不整合）：条件付き $f_{CAD}$ は0.0パーセンタイル（最も異常）と正しく検知します。一方、指標属性しか見ない無条件モデル $f_V$ は83.0パーセンタイルとむしろ正常寄りに評価してしまいます。$y$ の値自体は他のどこかの環境では普通に起こりうる値なので、環境を無視すると見逃されるのです。
つまり、環境と指標を同列に扱う手法は点Aで誤検知し、指標だけを見る手法は点Bで見逃す。両方の失敗モードを同時に回避できるのは、環境属性で条件づけた $f_{CAD}$ だけ——これが Song et al. (2007) が20年近く前に示した、条件付き異常検知の本質です。
最後に、6節の摂動プロトコルの考え方に沿って、$outliers$（環境だけが珍しい点）を実際に誤検知しないことも確認しておきます。
条件付き(f_CAD): perturbed recall=0.99 nonPerturbed誤検知率=0.01 outliers誤検知率=0.00
無条件(指標yのみ): perturbed recall=0.76 nonPerturbed誤検知率=0.23 outliers誤検知率=0.20
条件付きモデルは真の異常（perturbed）をほぼ完璧に検知しつつ、環境だけが珍しい正常点（outliers）を一切誤検知しません。無条件モデルは検知力そのものも劣るうえ、outliersの20%を誤って異常と判定してしまいます。論文のTable 1・Table 2が伝えたかった構図が、小さな合成データでも一貫して再現されました。
10. まとめ
Song, Wu, Jermaine, Ranka, “Conditional Anomaly Detection” (IEEE TKDE, 2007) を、原論文を精読しながらたどりました。
- 問題提起：異常検知が「候補の予算」を無駄にしないためには、環境的に珍しいだけの点を誤って報告しないことが重要。
- モデル：環境属性を条件とし、環境の混合ガウス $U$ ・指標の混合ガウス $V$ ・対応づけ $p(V_j\mid U_i)$ の3つ組で $f_{CAD}(y\mid\Theta,x) = \sum_i p(x\in U_i)\sum_j f_G(y\mid V_j)p(V_j\mid U_i)$ を構成する。
- 学習：EMアルゴリズムで、同時責任度 $b_{kij}$ を介して $U$・$V$・対応づけを更新する。実装上は3つの学習アルゴリズム（Direct-CAD／GMM-CAD-Full／GMM-CAD-Split）があり、実験ではGMM-CAD-Fullが最良。
- 評価：環境属性だけを操作した「摂動」による評価プロトコルで、13個の実データセットに対しCAD手法が単純GMM・$k$-NN・LOFを一貫して上回ることを示した。
- 実証：本記事の実装でも、環境と指標を同列に扱う結合GMMは環境が珍しいだけの点を誤検知し、指標だけを見るモデルは条件付き異常を見逃す一方、$f_{CAD}$ は両方を正しく判定することを数値で確認した。
この古典的な枠組みは、条件付きVAE・条件付き正規化フロー・Time-CADのような現代の文脈認識型異常検知手法すべての源流にあたります。次の記事では、より実務的な「運転条件の正規化」という、条件付きの考え方の別の実装を見ていきます。
主な参考文献
- X. Song, M. Wu, C. Jermaine, S. Ranka, “Conditional Anomaly Detection,” IEEE Transactions on Knowledge and Data Engineering, 19(5), 631–645, 2007.
- A. P. Dempster, N. M. Laird, D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” J. Royal Statistical Society B, 39(1), 1–38, 1977.
- E. M. Knorr, R. T. Ng, V. Tucakov, “Distance-Based Outliers: Algorithms and Applications,” VLDB Journal, 8(3-4), 237–253, 2000.（$k$-NN外れ値検知）
- M. M. Breunig, H.-P. Kriegel, R. T. Ng, J. Sander, “LOF: Identifying Density-Based Local Outliers,” ACM SIGMOD, 2000.
- C. M. Bishop, “Pattern Recognition and Machine Learning,” Springer, 2006.（混合モデルとEM）