Anomaly Detection

Normal Distribution

正規（ガウス）分布 Normal (Gaussian) distribution により、平均値を頂点として、データがどのように分散しているか（どのような確率分布を持っているか）を示すことができる。

${\scriptsize \text{$\mu = $ the mean of $x$}} \\ {\scriptsize \text{$\sigma^{2} = $ the variance of $x$ ($\sigma =$ standard deviation)}} \\ \begin{align} p(x; \mu, \sigma^{2}) & = \frac{1}{ \sqrt{ 2 \pi \sigma^{2} } } \exp \left( - \frac{ (x - \mu)^{2} }{ 2 \sigma^{2} } \right) \\ \end{align} \\$

μ は平均値
σ は標準偏差で、σ^2 で分散度合い（平均値からの距離）を示す。
確率密度関数 p をグラフにプロットすると、平均値を頂点として、左右対称のベル上の形状を描く。σ^2 値が大きいほど、頂点（平均値である確率）が下がり、勾配がなだらかになり、まんべんなく分布していることを示す。

Algorithm Evaluation

各パラメータ毎に正規分布から得られた確率値の積として、確率密度関数 p を定義する。

$\begin{align} \mu_{j} & = \frac{1}{m} \sum_{i = 1}^{m} x_{j}^{(i)} \\ \sigma_{j}^{2} & = \frac{1}{m} \sum_{i = 1}^{m} (x_{j}^{(i)} - \mu_{j})^{2} \\ p(x) & = p(x_1; \mu_1, \sigma_1^2) \cdot p(x_2; \mu_2, \sigma_2^2) \cdot \ldots \cdot p(x_n; \mu_n, \sigma_n^2) \\ & = \prod_{j = 1}^{n} p(x_{j}; \mu_{j}, \sigma_{j}^{2}) \\ & = \prod_{j = 1}^{n} \frac{1}{ \sqrt{ 2 \pi \sigma_{j}^{2} } } \exp \left( - \frac{ (x_{j} - \mu_{j})^{2} }{ 2 \sigma_{j}^{2} } \right) \\ \end{align} \\$

既知の正常データが 10,000 件、例外データが 20 件あるとすと、以下のように配分すると良い。確率分布を求める学習データは正常データのみである点に注意する。検証データとテストデータは同じ正常データでもよい。ただし例外データは分ける。

学習データ: 正常データのみの 6,000 件
検証データ: 正常データ 2,000 件 / 例外データ 10 件
テストデータ: 正常データ 2,000 件 / 例外データ 10 件

学習データから、確率密度関数 p（平均値 μ と分散度 σ^2）を決める。

次に、例外と見なす分岐点（確率）を決める。学習データから得られた確率密度関数より、検証データの確率値ベクトルを得る。確率値の最小値と最大値の間で、任意のステップ数で分岐点を取り出し、分岐点よりも確率値が小さければ、例外であるとみなす。

% pval: vector of probabilities
s = (max(pval) - min(pval)) / 1000;
for v in min(pval):s:max(pval)
    predictions = (pval < v); % vector of (0:normal|1:anomaly)
    ...
end

分岐点ごとに、精度 Precision と再現率 Recall を算出し F-score が最も良い分岐点を採用する。

${\scriptsize \text{$tp = $ true positive, $fp = $ false positive, $fn = $ false negative}} \\ {\scriptsize \text{$P = $ Precision, $R = $ Recall}} \\ P = \frac{tp}{tp + fp} \\ R = \frac{tp}{tp + fn} \\ F_{1} = 2 \frac{PR}{P + R} \\$

True Positive: 正常データを正常と判定した数（予測成功）
False Positive: 正常データを例外と判定した数（予測失敗）
False Negative: 例外データを正常と判定した数（予測失敗）

分岐点が決定したら、テストデータに対して例外を予測し、その精度を確かめる。

このように、確率分布から例外の分岐点を決める方法は Supervise learning と異なり、アルゴリズム（確率密度関数）の組み立てには正解値を必要としない。

現実に則した正常数と例外数（例：例外率 1 %）に分離されるような、入力パラメータ x を見つけることが重要になる。

正常時は比例しているパラメータが、異常時にそのバランスが崩れていることが多い。関連するパラメータを見つけ比率をとることで、例外時に突出した値になる。

不良 PC を見つけるとすると、一般的に CPU 利用率は処理数に比例するが、これらの値を個別にパラメータとして取っても、単に PC の貢献度を計るだけで、例外は検知できない。このような時は CPU 利用率とリクエスト処理数の比を取ることで、CPU が異常な使われ方をしていることをパラメータ化するとよい。

Multivariate Normal Distribution

正規分布では、パラメータ間の相関は持たない。例を上げると

CPU 利用率 10-90% を正常
秒間アクセス数 5-95 を正常

のように、範囲が各パラメータ空間内で閉じており、秒間アクセス数 5 で CPU 利用率 90% でも正常と判断するほかない。

例外とみなす分岐点は、各パラメータ内に限った上限と下限を決定することしかできず、二次元グラフにプロットした場合、分岐境界が描く円は、軸に対して平行 Axis-aligned になる。

秒間アクセス数 10 以下で CPU 利用率 80% 以上なら異常
CPU 利用率 20% 以下で、秒間アクセス数 50 以上なら異常

のように、CPU 利用率と秒間アクセス数が「比例していなれば例外」と判断したい時は、比率「CPU 利用率/アクセス数」を新しいパラメータとして追加する方法があるが、パラメータにどのような相関関係にあるかというのは、判断しづらい場合もある。

多変量正規分布 Multivariate normal distribution により、パラメータ間の相関関係も含めた多次元の正規分布を見ることができる。

${\scriptsize \text{$m = $ number of examples}} \\ {\scriptsize \text{$n = $ number of features}} \\ \begin{align} \mu & = \frac{1}{m} \sum_{i = 1}^{m} x^{(i)} \in \mathbb{R}^{n} \\ S & = \frac{1}{m} \sum_{i = 1}^{m} (x^{(i)} - \mu) (x^{(i)} - \mu)^{T} \in \mathbb{R}^{n \times n} \\ p(x; \mu, S) & = \frac{1}{ ( \sqrt{ 2 \pi } )^{m} \sqrt{ | S | } } \exp \left( - \frac{1}{2} (x - \mu)^{T} S^{-1} (x - \mu) \right) \\ \end{align}$

S は、入力ベクトル x と平均ベクトル μ の共分散行列 Covariance matrix
|S| は、行列式 Determinant = det(S)
S は逆行列をとっているため、可逆である必要がある。
S が正規分布の σ^2 の対角行列（対角成分以外が 0 で相関がない）ならば、正規分布と同じ値となる。
共分散から算出するため、パラメータ毎に分けた確率密度関数 p(j) を取る必要はない。

S が非可逆の特異行列 Singular matrix の場合には、以下の傾向がある。

サンプル数 m が、パラメータ数 n に対して少なすぎる。
重複しているパラメータ x1 = x2, x3 = x4 + x5, ... がある。

共分散行列から算出するため O(N^2) のコストとなる点が欠点になる。各パラメータに相関がないか、相関をパラメータ化できているならば、あえて使う必要はない。

vs. Normal Distribution

function p = normalDist(x, u, s)
    n = size(x, 1);
    ps = zeros(n, 1);
    for j = 1:n
      ps(j) = (1 / (sqrt(2 * pi * s(j)))) * exp(-1 * ((x(j) - u(j))^2) / (2 * s(j)));
    end
    p = prod(ps);
end

function p = multiNormalDist(x, u, S)
    n = size(x, 1);
    p = (1 / (((sqrt(2 * pi))^n) * sqrt(det(S)))) * exp(-1 * ((x - u)' * pinv(S) * (x - u)) / 2);
end

% CPU load, number of requests
X = [34 5; 56 10; 89 15; 90 20; 125 32; 68 18];

[m, n] = size(X);

% mean of X
u = mean(X)';
printf('mean: [%f, %f]\n', u);

% variance of X
s = (sum(((X - repmat(u', m, 1)) .^ 2)) ./ m)';
printf('sigma: [%f, %f]\n', s);

%
% Normal distribution:
%
normalDist([77; 16.6], u, s) % 6.5101e-04
normalDist([34; 5], u, s)    % 8.3081e-05
normalDist([125; 32], u, s)  % 3.1740e-05

% According to the results above, we assume that the probability
% in normal cases should be greater than 3.1740e-05.
%
% In this case, the CPU load is too high, although the number
% of requests is low. We expected that it would be labeled as
% an anomaly, but the probability with normal distribution is
% better than 3.1740e-05.
normalDist([100; 8], u, s) % 2.8090e-04

%
% Multivariate normal distribution:
%
v = X - repmat(u', m, 1);
S = (v' * v ./ m)';
multiNormalDist([77; 16.6], u, S) % 0.0020253
multiNormalDist([34; 5], u, S)    % 6.6206e-04
multiNormalDist([125; 32], u, S)  % 3.9053e-04

% The probability with multivariate normal distribution is
% sufficiently lower than the limit 3.9053e-04, so we can
% predict that this case must be an anomaly.
multiNormalDist([100; 8], u, S)   % 3.3423e-10

% If the covariance matrix "S" is a diagonal matrix of the
% sigma "s" of normal distribution, the probability will
% completely be the same as normal distribution.
multiNormalDist([77; 16.6], u, diag(s)) % 6.5101e-04
multiNormalDist([100; 8], u, diag(s))   % 2.8090e-04

vs. Supervised Learning