多変量統計的プロセス管理(MSPC)の数式と向きあってみる

はじめに

ひょんなことから多変量統計的プロセス管理 (MSPC: Multivariate Statistical Process Control) という異常検知の手法について調べる機会がありました。
MSPC の数式をまじめに説明している記事がネットにあまりなかったので記事にしてみました。

この記事には個人の意見や主観的な解釈が含まれます。間違っている点あればご指摘ください。

単変量統計的プロセス管理

MSPC の前に、MSPC の単変量版である単変量統計的プロセス管理(USPC: Univariate Statistical Process Control)について簡単に説明する。

USPC の代表例が健康診断らしい。
たとえば健康診断の血液検査では、白血球の数や γ-GTP などのそれぞれの項目に対して基準値が設定されているが、その基準値を超えた場合に異常と判断される。

USPC の問題点

引用：スモールデータ解析と機械学習 | Ohmsha

上記の図では、四角で囲まれた範囲は身長と体重の正常範囲を示している。
USPC ではこの四角の中にサンプルがあれば、正常であると判定される。

■ ◯ のサンプル
USPC ではいずれも正常と判定される

■ × のサンプル
USPC では異常とは判定されない！
× は肥満または痩せすぎの人たちを表しているが、この人たちは身長と体重の基準値を超えていないので正常体型であると判定される。

■ ☆ のサンプル
USPC では異常と判定される！
☆ の人は標準的な体型と比較すると大柄だが、身長と体重のあるべき関係からすると異常ではないので、異常値として判定されるのは不適切な場合がある。

つまり、USPC では複数の項目を複合的に加味した異常値の判定ができない！
(USPC は単変量の手法なので当たり前ではあるが)

多変量統計的プロセス管理 (MSPC) の概要

上記で異常と判定されなかった × のサンプルや、逆に異常と判定されてしまった ☆ のサンプルを意図通りに判定できるようにするための手法。

要するに、複数の項目の関係性(=相関係数)を考慮して異常値の判定を行う。

マハラノビス距離を使った異常検知と何が違う？

ここまで読んで、「それってマハラノビス距離でよくない？」と思った方はかなり鋭いです！

マハラノビス距離でも複数の項目の相関関係を考慮した異常検知はできるが、
MSPC は PCA(主成分分析)とマハラノビス距離を組み合わせた手法になっています。

T² 統計量と Q 統計量

PCA おさらい

PCA を行列で表現すると以下のようになる
(この式の導出は割愛します）

$$\mathbf{T} = \mathbf{X} \mathbf{V} \tag{1.1}$$

ここで、
$N$ : サンプル数
$P$ : 変数の数
$\mathbf{T}$ : 主成分得点行列 (サイズ : $N \times P$)
$\mathbf{X}$ : データ行列 (サイズ : $N \times P$)
$\mathbf{V}$ : $\mathbf{X}$の共分散行列の固有ベクトルを並べた行列 (サイズ : $P \times P$)
である。

主成分のうち上位$R$個のみを取り出した行列を $\mathbf{T}_R$, $\mathbf{V}_R$ とすると、(1.1)式は以下のようになる。

$$\mathbf{T}_R = \mathbf{X} \mathbf{V}_R \tag{1.2}$$

T² 統計量

(ちなみに T² はてぃーすくえあと呼ぶらしいです)

T² 統計量は以下で表される

$$T^2 = \sum_{r=1}^{R} \frac{t_r^{'2}}{\sigma_{t_r}^2} = \mathbf{t}_R^{'\top} \Sigma_R^{-2} \mathbf{t}_R^{'} \tag{2.1}$$

※ T² 統計量の T² は変数$T$の二乗ではなく、「T² 統計量」という統計量の名前である。

ここで、
$t_r^{'}$ : T² 統計量を計算したいデータに PCA を適用した際の、第$r$主成分得点 (スカラー)
$\sigma_{t_r}$: $\mathbf{T}_R$の列ベクトル(サイズ$N \times 1$)の標準偏差
$\mathbf{t}_R^{'}$: T² 統計量を計算したいデータに PCA を適用し、第$R$主成分までを取り出したベクトル (サイズ : $R \times 1$)
$\Sigma_R^2$ : $\mathbf{T}_R$の共分散行列 (サイズ : $R \times R$)

である。

マハラノビス距離の定義式は、

$$D = \sqrt{y^{\top} S^{-1} y}$$

$D$ : マハラノビス距離
$y$ : データ
$S$ : データの共分散行列

であるので、T² 統計量は、マハラノビス距離の一種と考えることができる。

(マハラノビス距離って何？という方：マハラノビス距離の意味を２次元の場合で理解する - 具体例で学ぶ数学)

もっと具体的にいうと、
T² 統計量はあるデータに PCA を行い、第$R$主成分までを取り出し他は捨てたデータで計算したマハラノビス距離と言える。

Q 統計量

Q 統計量は以下で表される

$$Q = \sum_{m=1}^{M} \left( x_m^{'} - \hat{x}_m^{'} \right)^2 = x^{'\top} \left( I - V_R V_R^{\top} \right) x^{'} \tag{3.1}$$

ここで、
$x^{'}$ : Q 統計量を計算したいデータ(サイズ : $P \times 1$)
$\hat{x}^{'}$ : $x^{'}$を PCA の第$R$主成分までで近似したデータ(サイズ : $P \times 1$)
$\mathbf{I}$ : 単位行列(サイズ : $P \times P$)
である。

$\left( x_m^{'} - \hat{x}_m^{'} \right)^2$の部分で分かるとおり、
Q 統計量はあるデータを表すベクトル$\mathbf{x}^{'}$と、それを PCA で近似したベクトル$\hat{\mathbf{x}}^{'}$の誤差の２乗和を表す指標である

ただ、$x^{'\top} \left( I - V_R V_R^{\top} \right) x^{'}$の部分がよくわからないと思うので、
前準備として上記の PCA の式展開をもう少し行う。

(1.2)式より、

$$\mathbf{X} = \mathbf{T} \mathbf{V}^{\top} \tag{3.2}$$

PCA において、第$r$主成分得点まで取り出した $\mathbf{T}_R$, $\mathbf{V}_R$で$\mathbf{X}$を近似すると、(3.2)式は以下となる

$$\mathbf{X} \approx \hat{\mathbf{X}} = \mathbf{T}_R \mathbf{V}_R^{\top} \tag{3.3}$$

また、(1.2)式と(3.3)式より、

$$\hat{\mathbf{X}} = \mathbf{X} \mathbf{V}_R \mathbf{V}_R^{\top} \tag{3.4}$$

続いて、$\hat{\mathbf{X}}$で近似しきれなかった部分を $E$ とすると、

$$E = \mathbf{X} - \hat{\mathbf{X}} = \mathbf{X} - \mathbf{X} \mathbf{V}_R \mathbf{V}_R^{\top} \tag{3.5} \\ = \mathbf{X} \left( \mathbf{I} - \mathbf{V}_R \mathbf{V}_R^{\top} \right)$$

ここで、$\mathbf{V}_R \mathbf{V}_R^{\top}$について考えてみる。

実は、$\mathbf{V}_R \mathbf{V}_R^{\top}$は直交射影行列になっている。

(直交射影行列って何？という方：射影行列のイメージと楽しい公式 | 高校数学の美しい物語)

直交射影行列の証明には、$\mathbf{A}\mathbf{A} = \mathbf{A}$となること、 $\mathbf{A}^{\top} = \mathbf{A}$となることをそれぞれ示せば良い。

$$\begin{aligned} \mathbf{V}_R \mathbf{V}_R^{\top} \mathbf{V}_R \mathbf{V}_R^{\top} &= \mathbf{V}_R \mathbf{V}_R^{\top} \\ \left( \mathbf{V}_R \mathbf{V}_R^{\top} \right)^{\top} &= \mathbf{V}_R \mathbf{V}_R^{\top} \\ \end{aligned}$$

($\mathbf{V}_R$は固有ベクトルを並べた行列で、固有ベクトル同士は直交するので$\mathbf{V}_R^{\top} \mathbf{V}_R = \mathbf{I}$です)

両方とも成り立つので、確かに$\mathbf{V}_R \mathbf{V}_R^{\top}$は直交射影行列である。

さらに、直交射影行列のもつ性質として
ある直交射影行列$A$において$\mathbf{I} - \mathbf{A}$も直交射影行列であり、かつ$\mathbf{A}$の直交補空間へ射影する行列になる。

(直交補空間って何？という方：直交補空間の性質 | 高校数学の美しい物語)

つまり、以下が成り立つ。

$$\mathrm{Im}(\mathbf{V}_R \mathbf{V}_R^{\top}) = \mathrm{Im}(\mathbf{I} - \mathbf{V}_R \mathbf{V}_R^{\top})^{\perp}$$

補足：$\perp$は直交補空間を表し、$\mathrm{Im}(\cdot)$は行列によって射影される線形空間を表す。

よって、(3.4), (3.5)式は以下のように解釈できる。

■ (3.4)式：

$$\hat{\mathbf{X}} = \mathbf{X} \mathbf{V}_R \mathbf{V}_R^{\top}$$

$\mathbf{X}$を$\hat{\mathbf{X}}$に近似することは、$\mathbf{X}$を構成する行ベクトル$\mathbf{x}$を$\mathbf{V}_R \mathbf{V}_R^{\top}$で正射影することに相当する。

さらに言うと、
(1.2)式より$\mathbf{V}_R$はデータ行列$\mathbf{X}$を主成分得点$\mathbf{T}_R$に変換する行列なので、
$\mathbf{X} \mathbf{V}_R \mathbf{V}_R^{\top}$は一回$\mathbf{X}$を主成分の軸で構成される空間に写像し、それをまた元の空間に写像する変換とみなせる。
（ただし$\mathbf{V}_R$は上位$R$個の主成分のみの行列のため、元通りの$\mathbf{X}$にはならない）

■ (3.5)式：

$$E = \mathbf{X} - \hat{\mathbf{X}} = \mathbf{X} - \mathbf{X} \mathbf{V}_R \mathbf{V}_R^{\top} \\ = \mathbf{X} \left( \mathbf{I} - \mathbf{V}_R \mathbf{V}_R^{\top} \right)$$

$\mathbf{X}$を$\hat{\mathbf{X}}$に近似した際の誤差$E$は、$\mathbf{X}$を構成する行ベクトル$\mathbf{x}$を$(\mathbf{I} - \mathbf{V}_R \mathbf{V}_R^{\top})$で正射影することに相当する。

前準備が終わったので、再び Q 統計量の式に戻る。

式(3.1)を$\left( I - V_R V_R^{\top} \right)$が射影行列であることを利用し変形すると、以下のようになる。

$$Q = \sum_{m=1}^{M} \left( x_m^{'} - \hat{x}_m^{'} \right)^2 = x^{'\top} \left( I - V_R V_R^{\top} \right) x^{'} \\ = x^{'\top} \left( I - V_R V_R^{\top} \right) \left( I - V_R V_R^{\top} \right) x^{'} \\ = \left( \left( I - V_R V_R^{\top} \right) x^{'} \right)^{\top} \left( I - V_R V_R^{\top} \right) x^{'}$$

つまり、Q 統計量は$\left( I - V_R V_R^{\top} \right) x^{'}$、つまり誤差$E$の二乗ノルムとみなすことができる。

T² 統計量と Q 統計量の気持ちになって考える

T² 統計量は前述のとおり、あるデータに PCA を行い第$R$主成分までを取り出し他は捨てたデータで計算したマハラノビス距離である。
マハラノビス距離は単体でも異常検知手法として使われることがあるため、T² 統計量を異常を検知する指標として使うのはとても自然であると言える。

ただし、T² 統計量は第$R$主成分までのみを使うということで、ここには暗黙的に「第$R$主成分まででデータの代表的な特徴を捉えられる」という仮定が含まれている。

なので上記の仮定が崩れる、つまりデータの質が変わり PCA による近似がうまくいかなくなった場合は、それ自体も異常として検知したい。
これを捉えるために有効なのが Q 統計量で、Q 統計量の実態としては PCA による近似誤差の二乗和である。

「PCA の第$R$主成分までのみを使う」というのは「データ行列$\mathbf{X}$を$\hat{\mathbf{X}}$とみなす」のと数式上等価であるのが、MSCP の数式を理解する上でのキモに感じました（個人的な解釈）
そのため、PCA の近似誤差は$\mathbf{X}$と$\hat{\mathbf{X}}$の差で計算できる。

結局、T² 統計量と Q 統計量どう使えばいいのか？

引用： 研究とベンチャーと数学と（その３）｜ Dr. Kano

MSPC では、まずデータを PCA にかけて、代表的な主成分とそれ以外の雑音の２つの成分にデータを分割する。

なぜ PCA を使うかというと、データのもつ「異常」を以下２つの種類に分けて検知したいため。

1. 監視対象のデータの代表的な相関関係に反している異常
2. 監視対象のデータの代表的ではない相関関係 (一部のサンプルのみが関係していたり、ただのノイズによって生まれた相関) に反している異常

MSCP では、前者を$T^2$統計で、後者を$Q$統計量を使って監視する。

個人的にはこれを読んだとき「$Q$統計量ってほんとにいる？」と思ったので
$Q$統計量も監視した方がいいと言える具体例がほしいなと思ったが、ネットにあまり落ちていないので 最終手段として ChatGPT に聞いてみた。

ChatGPT への質問:
$T^2$統計量では異常が検出されず、$Q$統計量では異常が検出される例を具体的なシーンを例にあげて３つほど紹介して

回答(ChatGPT):
T² 統計量では異常が検出されず、Q 統計量では異常が検出されるシーンは、主にプロセスが複雑で、全体的な平均の変化は見られないが、プロセスの特定の側面で異常が生じるケースに関連しています。以下に具体例を 3 つ挙げます。

1. 製造プロセスでの特定機器の故障

• シーン: 複数の製造ラインが同時に稼働しており、それぞれが製品の異なる部分を担当しています。プロセス全体の平均的なパフォーマンスには大きな変化はありませんが、ある一つの機器が不具合を起こし、その機器に関係する特定の工程だけが異常を示します。
• 結果: T² 統計量はプロセス全体をモニタリングしているため、平均的な変動がない場合、異常を検出しません。しかし、Q 統計量は PCA モデルで説明できない変動（特定の機器の故障による局所的な異常）を捉えるため、異常が検出されます。

2. 品質管理での特定成分の異常

• シーン: 食品製造において、複数の成分がミックスされる過程があります。全体的な成分バランスは正常に見えるものの、ある特定の成分の投入量がわずかにずれ、最終製品の品質に影響を与える可能性があります。このずれは他の成分のバランスで相殺され、全体的な成分バランスは保たれています。
• 結果: T² 統計量は全体のバランスに基づいて異常を検出しないが、Q 統計量は PCA モデルが捉えられなかったその特定成分の異常な挙動を検出します。

3. センサー故障による異常データ

• シーン: 工場の生産ラインで複数のセンサーが温度や圧力などを監視しています。全体的なプロセスのパフォーマンスは正常であるが、特定のセンサーが故障し、不正確なデータを報告しています。この故障により、一部のデータが通常のプロセス挙動から外れることになります。
• 結果: T² 統計量は全体的なデータの平均値が大きく変わらない限り異常を検出しませんが、Q 統計量は PCA モデルによって説明されないセンサーの異常なデータを検出します。

なんとなく、、雰囲気は理解できたような、、、気がする、、、

まとめ

異常検知手法の１つである MSPC について、数式メインで説明してみました。
MSPC は 1979 年からある古い手法ではありますが、変数の数が多いデータや PCA がうまくハマるデータに対しては有効な手法であると感じました。

あと数式を眺めると線形代数の良い復習になってよかったです。