入門 機械学習による異常検知―Rによる実践ガイド： ノート3

本読み 機械学習

読んでいる本（出典）： 入門 機械学習による異常検知―Rによる実践ガイド | 井手 剛 |本 | 通販 | Amazon

前回： ノート2 ／ 次回： まだ
目次： 入門 機械学習による異常検知―Rによる実践ガイド

読んだページ： 44～52ページ
以下、メモと雑談。

前回までのあらすじ

  • データが多次元正規分布にしたがう場合、 a(x') = (x' - \hat{\mu}) ^{\rm T} \hat{\Sigma} ^{-1} (x' - \hat{\mu}) で異常度を定義できる。
    観測データが i.i.d. に多次元正規分布にしたがうなら、この異常度は F 分布にしたがう。
    • 例．（ノート2と同じ例）2017年1月1日～18日の東京の最低／最高気温は、1.63±1.99℃／10.47±2.77℃だった。
      気象庁｜過去の気象データ検索
      このうち、2017年1月8日の東京の最低気温は1.6℃、最高気温は6.0℃だった。2017年1月8日は異常な日だったのだろうか。但し、ここで恣意的に、「14%も起こらないこと」を異常の閾値にする。
      → A. 最低気温と最高気温を1変数ずつみると異常ではないが、2変数の組合せは異常となる（以下）。
      • 2017年1月8日の最低気温だけみると1.6℃で、これはこの18日間の平均値1.63℃にかなり近いので異常値ではないだろう（証明略）。
      • 2017年1月8日の最高気温だけみると6.0℃とこの18日間で2番目に低いが、最高気温が独立同一正規分布にしたがうなら、「14%も起こらないこと」ではない（異常度 2.331211 < F分布の上側14%点 2.396771）。
        x1 <- c(2.0, 3.8, 3.5, 3.6, 3.7, 1.5, 0.1, 1.6, 3.8, 3.5, 3.9, 0.7, 1.4, -1.3, -2.3, -2.0, 0.7, 1.1)
        x2 <- c(13.8, 13.3, 13.7, 14.0, 10.4, 8.8, 8.7, 6.0, 11.1, 12.7, 11.0, 12.1, 12.7, 6.3, 4.7, 8.0, 10.9, 10.3)
        > a <- (18-1)/(18+1) * (x2[8]-mean(x2))^2 / (var(x2)*(18-1)/18)
        > a
        [1] 2.331211
        > F_0.14 <- qf(df1=1, df2=18-1, 0.86)
        > F_0.14
        [1] 2.396771
      • 2017年1月8日の最低／最高気温の組合せだと、最低／最高気温が独立同一2次元正規分布にしたがうなら、「14%も起こらないこと」といえる（F分布の上側14%点 2.228783 < 異常度 2.248451）。
        > mu <- colMeans(cbind(x1, x2))
        > xc <- cbind(x1, x2) - matrix(1, 18, 1) %*% mu
        > sigma <- t(xc) %*% xc / 18
        > a <- rowSums( (xc %*% solve(sigma) ) * xc)
        > a <- a * (18 - 2) / ((18 + 1) * 2)
        > a[8]
        [1] 2.248451
        > F2_0.14 <- qf(df1=2, df2=18-2, 0.86)
        > F2_0.14
        [1] 2.228783

  • なぜ異常度が F 分布にしたがうかという以前に、i.i.d. に M 次元正規分布にしたがう観測データからの、平均ベクトルと共分散行列の最尤推定値の導出の復習…。
    • 対数尤度は \displaystyle L(\mu, \Sigma | x_{1:N}) = - \frac{MN}{2} \ln (2\pi) -\frac{N}{2} \ln |\Sigma| - \frac{1}{2} \sum_{n=1}^{N}(x_n - \mu)^{\rm T} \Sigma^{-1}(x_n - \mu) 。これを \mu と \Sigma について偏微分しないといけない。
      • \mu について偏微分するには、行列の積のトレースの微分公式 \displaystyle \frac{\partial}{\partial A} {\rm Tr} (ABA^{\rm T}) = (B + B^{\rm T})A をつかうと \displaystyle \frac{\partial L }{\partial \mu} = - \frac{1}{2} \sum_{n=1}^{N} \frac{\partial}{\partial \mu} {\rm Tr} \Bigl( (x_n - \mu)^{\rm T} \Sigma^{-1}(x_n - \mu) \Bigr) = - \Sigma^{-1} \sum_{n=1}^{N} (x_n - \mu)
        （∵ 対称行列 \Sigma は逆行列もまた対称なので、 \Sigma^{-1} + (\Sigma^{-1})^{\rm T} = 2\Sigma^{-1}）
      • \Sigma について偏微分するには、正則な正方行列に関する微分公式 \displaystyle \frac{\partial \ln |A|}{\partial A} = (A^{-1})^{\rm T} と、
        ベクトル a と対称行列 B について成り立つ a^{\rm T} B a = \sum_i \sum_j b_{i,j} a_i a_j = {\rm Tr}(B a a^{\rm T}) と、
        行列の積のトレースの微分公式 \displaystyle \frac{\partial}{\partial B} {\rm Tr} (BA) = A^{\rm T} をつかうと
        \displaystyle \frac{\partial L}{\partial (\Sigma^{-1})} = \frac{N}{2} \frac{\partial}{\partial (\Sigma^{-1})} \ln |\Sigma ^{-1}| - \frac{1}{2} \sum_{n=1}^{N} \frac{\partial}{\partial (\Sigma^{-1})} (x_n - \mu)^{\rm T} \Sigma^{-1}(x_n - \mu)
        　　 　　 \displaystyle \, = \frac{N}{2} \Sigma - \frac{1}{2} \sum_{n=1}^{N} (x_n - \mu) (x_n - \mu)^{\rm T}
  • それで、なぜ異常度 a(x') = (x' - \hat{\mu}) ^{\rm T} \hat{\Sigma} ^{-1} (x' - \hat{\mu}) が F 分布にしたがうのか（44～48ページ）。
    • 証明の仕方は技巧的（47ページ）で、 \displaystyle a(x') = \frac{(x' - \hat{\mu}) ^{\rm T} \Sigma ^{-1} (x' - \hat{\mu})}{(x' - \hat{\mu}) ^{\rm T} \Sigma ^{-1} (x' - \hat{\mu}) / (x' - \hat{\mu}) ^{\rm T} \hat{\Sigma} ^{-1} (x' - \hat{\mu})} とすると分子も分母もカイ2乗分布にしたがうので、全体はF分布にしたがう、という仕組み。
    • 分母がカイ2乗分布にしたがう説明： 「分母は。もともと自由度 N-1 のウィシャート分布にしたがう M \times M 行列を、ブロック分割行列の逆行列の公式（付録の定理 A.4 参照）を使って 1 \times 1 行列に縮約したものに対応している（48ページ）」とあるけど、定理 A.4 を参照してもここの意味がよくわからない…。なぜ逆行列の公式で縮約ができるのか…。線型代数は大学で習ったけどもう覚えていなかった…。
    • ウィキペディアには証明はなかった…原論文へのリンクは有。

今回のお話

  • 多次元のホテリング理論で「最高／最低気温の組合せが異常」というのは検知できるが、もっといえば、最高気温と最低気温のどちらが／あるいは両方が異常だったのかも知りたい。→ 正常であることが期待されるデータ群 \mathcal{D} から、全データの1変数あたりのマハラノビス距離を計算し、ここまでが正常という閾値を適当に決め、異常なデータの1変数あたりのマハラノビス距離をこの閾値と比較する。