Hotelling の T^2 で異常値にラベルを付ける

最終更新 2025-11-03 読了時間 3 分

まとめ

Hotelling の T^2 統計量は、各観測値が平均からどれだけ離れているかを分散で割って測るマハラノビス距離の 2 乗です。データが多変量正規分布に従うなら T^2 は自由度 $d$ の $\chi^2$ 分布に従うので、scipy.stats.chi2.ppf(1 - alpha, df=d) で「正常なら 95% がこの値以下に収まる」という閾値が 1 行で決まります。このページでは正規乱数 500 点に 10 個の異常値を混ぜた時系列を作り、T^2 が閾値 3.84 を超えた点に異常ラベルを付けます。読み終えると、alpha を 0.05 から 0.01 に下げると検出される点がどう減るか、1 次元の T^2 が z スコアの 2 乗と同じものだと説明できるようになります。

直感 #

工場のセンサーが毎秒 1 つの値を返しているとします。ほとんどの値は平均 0 のまわりに標準偏差 1 でばらつき、たまに 6 前後の値が飛び込みます。「平均から標準偏差の何倍離れたか」を 2 乗した値が T^2 で、1 次元では z スコアの 2 乗そのものです。正常な点の T^2 は $\chi^2_1$ 分布に従うので、「正常な点のうち 5% しか超えない値」= 3.84 を線として引けば、それを超えた点に機械的に異常ラベルを付けられます。多変量になっても考え方は同じで、分散で割る代わりに共分散行列 $\Sigma^{-1}$ をはさみ、$\chi^2$ の自由度を次元数 $d$ に変えるだけです。

式で書くと何を測っているか #

特徴量ベクトル $x \in \mathbb{R}^d$、平均 $\mu$、共分散 $\Sigma$ に対し、

$$ T^2 = (x - \mu)^\top \Sigma^{-1} (x - \mu) $$

を計算します。$(x - \mu)$ が「平均からのずれ」、$\Sigma^{-1}$ が「ばらつきの大きい方向ほど小さく数える重み」です。多変量正規を仮定すると $T^2$ は自由度 $d$ の $\chi^2_d$ に従うため、検出閾値 $\tau$ は

$$ \tau = \chi^2_d(1 - \alpha) $$

で与えられます。$T^2 > \tau$ となる観測は異常と判定できます。1 次元では $\Sigma$ が分散 $\sigma^2$ 1 個になるので、T^2 は z スコア $(x-\mu)/\sigma$ の 2 乗と一致し、$\alpha = 0.05$ の閾値は $\chi^2_1(0.95) = 3.84$、つまり「z スコアの絶対値が 1.96 を超えたら異常」と同じ判定になります。

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：少ない回数で孤立する点ほど異常。外れ値は 1〜5 回で孤立し、密集した点は 8 回切っても 5 点が残る。T^2 も「密集した点から遠い点を異常とみなす」点は同じで、違いは距離を $\chi^2$ 分布の閾値 3.84 と比べて判定するところです。

サンプルデータを作成する #

次のコードは平均 0・標準偏差 1 の正規乱数 500 点を作り、ランダムに選んだ 10 点に平均 6 のずれを足して異常値にします。正解を知っているデータなので、後で検出結果と赤い × が重なるかを確認できます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(100)
n_samples = 500
time = np.arange(n_samples)
series = rng.normal(loc=0.0, scale=1.0, size=n_samples)

anomaly_index = rng.choice(n_samples, size=10, replace=False)
series[anomaly_index] += rng.normal(loc=6.0, scale=1.0, size=anomaly_index.size)

plt.figure(figsize=(10, 3))
plt.plot(time, series, ".", alpha=0.8, label="観測値")
plt.plot(anomaly_index, series[anomaly_index], "x", ms=10, c="red", label="挿入した異常値")
plt.title("合成した時系列データ")
plt.xlabel("時刻")
plt.ylabel("値")
plt.legend()
plt.tight_layout()
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

4 行目乱数の種を固定する
default_rng(100) で種を固定しているので、異常値の位置も検出結果も実行ごとに変わりません。種を外すと図の赤い × の位置が毎回変わり、本文の「10 個中いくつ検出できたか」の説明と合わなくなります。
7 行目正常な点は平均 0・標準偏差 1
loc=0.0, scale=1.0 の正規乱数が正常な観測値です。この分布を仮定しているからこそ、後で T^2 を $\chi^2$ 分布と比べられます。実データで裾が重い場合は、この仮定が崩れて誤検知が増えます。
9-10 行目10 点だけ平均 6 のずれを足す
rng.choice(..., size=10, replace=False) で重複なく 10 個の時刻を選び、そこに平均 6・標準偏差 1 のずれを足しています。正常値の標準偏差が 1 なので、ずれ 6 は z スコアでおよそ 6、T^2 ではおよそ 36 に相当し、閾値 3.84 を大きく超える「見つけやすい異常」です。loc を 2 程度に下げると正常値の裾に埋もれ、取り逃がしが出始めます。
13-14 行目正解の位置を赤い × で重ねる
観測値を点、挿入した異常の時刻を赤い × で重ねています。次の図で「検出された異常値」の赤丸がこの × と重なれば、T^2 の判定が正しく働いたことになります。

サンプルデータを作成するの図

T^2 を計算してラベル付けする #

次のコードは系列の平均と分散から各点の T^2 を求め、chi2.ppf(0.95, df=1) = 3.84 を閾値にして超えた時刻を取り出します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
from scipy.stats import chi2

def detect_hotelling_t2(series: np.ndarray, alpha: float = 0.05) -> tuple[np.ndarray, np.ndarray, float]:
    series = np.asarray(series)
    mean = series.mean()
    var = series.var(ddof=1)
    t2_scores = ((series - mean) ** 2) / var
    threshold = chi2.ppf(1 - alpha, df=1)
    detections = np.where(t2_scores > threshold)[0]
    return detections, t2_scores, threshold

detected_index, t2_scores, threshold = detect_hotelling_t2(series, alpha=0.05)

plt.figure(figsize=(10, 3))
plt.plot(time, t2_scores, label="T^2 スコア")
plt.axhline(threshold, color="crimson", linestyle="--", label=f"閾値 (α = 0.05)")
plt.scatter(detected_index, t2_scores[detected_index], c="red", marker="o", s=60, label="検出された異常値")
plt.xlabel("時刻")
plt.ylabel("T^2")
plt.title("Hotelling の T^2 スコアと検出結果")
plt.legend()
plt.tight_layout()
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

5-6 行目平均と不偏分散を推定する
mean() と var(ddof=1) が式の $\mu$ と $\Sigma$ の 1 次元版です。異常値を含めたまま推定しているので、異常が多いデータでは分散が膨らんで閾値が甘くなります。その場合は中央値と MAD で置き換えるか、異常を除いた期間で推定し直します。
7 行目T^2 は z スコアの 2 乗
((series - mean) ** 2) / var は $(x-\mu)^2/\sigma^2$、つまり z スコアの 2 乗です。多変量に広げるときはこの行を (x - mu) @ inv(Sigma) @ (x - mu) に置き換えます。2 乗しているので、平均より上のずれも下のずれも同じ向きの大きな値になります。
8 行目閾値を χ² 分布から取る
chi2.ppf(1 - alpha, df=1) は「自由度 1 の $\chi^2$ 分布で上側 5% になる値」で 3.84 を返します。alpha=0.01 にすると 6.63 に上がり、検出が厳しくなって誤報が減る代わりに小さな異常を見逃します。多変量なら df を次元数に変えます。
9 行目閾値超えの時刻を取り出す
np.where(t2_scores > threshold)[0] が異常ラベルそのものです。この配列が空なら異常なし、長さが 10 前後なら挿入した異常をほぼ拾えています。正常値の 5% も確率的に閾値を超えるので、500 点中 25 点程度の誤検知が混ざるのは仕様どおりです。
15-17 行目スコアと閾値線を同じ図に描く
T^2 を折れ線、閾値を赤い破線、検出点を赤丸で重ねています。赤丸が前の図の × と同じ時刻にあれば検出成功、× のない場所の赤丸は誤検知です。alpha を変えて破線の高さを動かすと、赤丸の数が変わる様子を比べられます。

T^2 を計算してラベル付けするの図

実行結果の読み方 #

挿入した 10 点の T^2 はおよそ 25〜50 に達し、閾値 3.84 をはるかに超えるので赤丸で拾われます。一方、正常な 490 点のうちおよそ 5%、つまり 25 点前後も閾値をわずかに超えて赤丸になります。これは $\alpha = 0.05$ が「正常でも 5% は異常扱いにする」という設定だからで、誤報を減らしたいなら alpha=0.01（閾値 6.63）に下げます。

使いどころと落とし穴 #

次元が高い場合は共分散行列 $\Sigma$ の逆行列を直接計算せず、scipy.linalg.solve で連立方程式として解くと数値的に安定します。閾値を決める $\alpha$ は検出したい異常率に合わせて調整し、アラートを減らしたい監視用途では 0.01 など小さい値から試します。正規性が崩れる右に裾を引いたデータでは、先に Box-Cox や Yeo-Johnson 変換で対称化するか、RobustScaler のように中央値ベースの推定に替えると誤検知が減ります。オンライン監視では平均と共分散を一定期間ごとに再計算して、データドリフトに追随させてください。分布を仮定したくない場合や 2 次元以上で形が複雑なデータでは、次のページの Isolation Forest が代替になります。