chap3_Introduction_to_Anomaly_Detection_using_Machine_Learning

3161 days ago by takepwave

Hiroshi TAKEMOTO (take@pwv.co.jp)

入門機械学習による異常検出

井出　剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。

3章　非正規データからの異常検知

井出本の基本は、データに対するモデルを使って負の対数尤度を求め、それを異常検出関数として使うことです。

分布が左右対称でない場合

Davisの体重（weight）のヒストグラムを見ると、分布が左右対称ではないことに気づきます。 井出本では、ガンマ分布を使ってこの非対称な分布の異常度を検出しています。

準備

いつものように必要なパッケージを読込ます。

# RとPandasで必要なユーティリティ # Rの必要なライブラリ r('library(ggplot2)') r('library(jsonlite)') # RUtilにRとPandasのデータフレームを相互に変換する関数を追加+GgSaveを追加（2015/07/12） load(DATA + 'RUtil.py') 
       
# 2.2 carパッケージのDavisを使って例題を試す r('library(car)') r('library(MASS)') 
        
 [1] "MASS"      "car"       "jsonlite"  "ggplot2"   "stats"    
"graphics"  "grDevices" "utils"    
 [9] "datasets"  "methods"   "base"     
 [1] "MASS"      "car"       "jsonlite"  "ggplot2"   "stats"     "graphics"  "grDevices" "utils"    
 [9] "datasets"  "methods"   "base"

体重のヒストグラム

Davisの体重のヒストグラムを再度プロットしてみます。 1個だけ160のところにありますが、2章ではこれが体重と身長を入れ間違えたものと結論づけています。

# Rのデータをpandaのデータフレームに変換する davis = RDf2PandaDf("Davis") # weightの分布をみる # 1個だけ、160以上の値を示している ggplot(davis, aes(x='weight')) + geom_histogram(binwidth=5.0,color="grey") 
        
<ggplot: (8738330208357)>
<ggplot: (8738330208357)>
GgSave('fig_2.1.png', fac=0.8) 
        
Saving 11.0 x 8.0 in image.

Saving 11.0 x 8.0 in image.

ガンマ分布の対数尤度

2章と同様にガンマ分布の対数尤度を求めてみます。

井出本では、ガンマ分布を以下の様に定義しています。 $$ \mathcal{G} (x | k, s) = \frac{1}{s \Gamma(k)}(\frac{x}{s})^{k-1} exp(- \frac{x}{s}) $$

これから対数尤度Lは、以下の様になります。 $$ L(k, s | \mathcal{D}) = \sum_{n=1}^N \left [ -ln\{s \Gamma(k) \} + (k-1) ln\frac{x^{(n)}}{s} - \frac{x^{(n)}}{s} \right ] $$

この式をパラメータkとsで微分し、０と等しいとします。

$-ln\{s \Gamma(k) \}$が$-ln(s) - ln(\Gamma(k))$となりますから、sの偏微分は以下の様になります。 $$ \begin{eqnarray} \frac{\partial L}{\partial s} & = & \sum_{n=1}^N \left [ -\frac{1}{s} + (k-1)\frac{\partial}{\partial s} ln \frac{x^{(n)}}{s} + \frac{x^{(n)}}{s^2} \right ] \\ & = & \sum_{n=1}^N \left [ -\frac{1}{s} - \frac{(k-1)}{s} + \frac{x^{(n)}}{s^2} \right ] \\ & = & \sum_{n=1}^N \left [ -\frac{k}{s} + \frac{x^{(n)}}{s^2} \right ] \end{eqnarray} $$ この式が0となることから、 $$ \hat{s} = \frac{1}{\hat{k} N} \sum_{n=1}^N x^{(n)} = \frac{\hat{\mu}}{\hat{k}} $$

同様にkの偏微分は以下の様になります。 $$ \frac{\partial L}{\partial k} = \sum_{n=1}^N \left [ -\frac{1}{\Gamma(k)}\frac{\partial \Gamma(k)}{\partial k} + ln \frac{x^{(n)}}{s}\right ] $$ しかし、kがガンマ関数が入っているので、簡単には$\hat{k}$は求まりません。

モーメント法を使ったk, sの推定

井出本では、最尤推定でk, sを求める代わりにモーメント法を使ってkとsを推定しています。 1次のモーメントは以下の様に計算されます。$\Gamma(k+1) = k \Gamma(k)$を関係を使っています。(以下の式は未フォロー) $$ \lt x \gt = \int_0^{\infty} dx \, x \, \mathcal{G}(x | k,s) = s\frac{\Gamma(k+1)}{\Gamma(k)} \int_0^{\infty} dx \, x \, \mathcal{G}(x | k+1,s) = ks $$ 同様に２次のモーメントは、以下の様になります。 $$ \lt x^2 \gt = k(k+1)s^2 $$

これとデータから求められる１次の２次モーメントの関係 $$ \lt x \gt = \frac{1}{N} \sum_{n=1}^N x^{(n)} $$ $$ \lt x^2 \gt = \frac{1}{N} \sum_{n=1}^N {x^{(n)}}^2 $$

この関係から $$ \hat{k}_{mo} = \frac{\hat{\mu}^2}{\hat{\sigma}^2}, \hat{s}_{mo} = \frac{\hat{\sigma}^2}{\hat{\mu}} $$

davisの体重を使って計算してみる

davisのweight(体重）でモーメント法を使ってk, sを求めてみます。

# ガンマ分布の異常検出 N = len(davis.weight) mu = davis.mean()['weight'] si = davis.std()['weight']*(N-1)/N kmo = (mu/si)^2 smo = si^2/mu print kmo, smo 
        
19.1928236809 3.42836474164
19.1928236809 3.42836474164

Rのfitdistrを使ってガンマ分布を推定

Rのfitdistrを使ってdavisのweight(体重）のガンマ分布を推定してみます。 shapeがkの値、rateが1/sの値として返されます。

モーメント法とfitdistrの結果は、そこそこ近い結果となっています。

r('ml <- fitdistr(Davis$weight, "gamma")') 
        
     shape         rate   
  22.4854793    0.3417247 
 ( 2.2317648) ( 0.0342978)
     shape         rate   
  22.4854793    0.3417247 
 ( 2.2317648) ( 0.0342978)
kml = r('ml$estimate["shape"]') sml = 1/r('ml$estimate["rate"]') print kml print sml 
        
   shape 
22.48548 
    rate 
2.926333 
   shape 
22.48548 
    rate 
2.926333

scipyのstatsを使ってガンマ分布を推定

このノートの特徴は、Sageとpythonを使って井出本と同じ結果を求めることにあります。

scipyのstatsのgamma.fitを使って、kとsを求めてみます。 gamma.fitでは、floc=0を指定しないとfitdistrと同じ結果になりません。

# SciPyのstatsを使ってガンマ分布の最尤推定値を求める、floc=0を指定するとfitdistrと同じ結果になる from scipy import stats fit_alpha, fit_loc, fit_beta = stats.gamma.fit(np.array(davis.weight), floc=0) print(fit_alpha, fit_loc, fit_beta) 
        
(22.485303592007426, 0, 2.9263558630975797)
(22.485303592007426, 0, 2.9263558630975797)

結果の可視化

今回バージョンアップしたSage6.7からヒストグラムがサポートされました。

通常のヒストグラムと正規化（normed=True）したヒストグラムが表示でき、 他のプロットデータと重ね合わせることができます。

# Pandasのデータフレームからweightを取り出す weight = list(davis.weight.values) 
       
# Sage6.7で導入されたhistogramを使用、binsのデフォルトは10 histogram(weight, figsize=5, bins=24) 
       
# 正規化したグラフも可能 hg = histogram(weight, figsize=5, bins=24, normed=True) hg.show() 
       
def _gamma(x, k, s): return 1/(s*gamma(k))*(x/s)^(k-1)*exp(-x/s) 
       
gp = plot(_gamma(x, 22.48548, 2.926333), [x, 0, 180], rgbcolor='red') (gp+hg).show() 
       

2章の正規分布との重ね合わせ

以下の様に正規分布関数_gaussを定義して、ヒストグラム、正規分布、ガンマ分布の推定結果を 重ね合わせてみます。

def _gauss(x, mu, sig2): return 1/sqrt(2*pi*sig2)*exp(-1/(2*sig2)*(x - mu)^2) 
       
# 正規分布として推定した結果（2章の結果を利用） sig2 = 226.7 mu = 65.8 ga = plot(_gauss(x, mu, sig2), [x, 0, 180], rgbcolor='green') (hg+gp+ga).show(figsize=5) 
       

stats.gamma.fitの精度

最初、flocを指定しないでフィッティングをしたため、fitdistrと結果が異なり 悩んでいました。しかしfloc=0とすることで、fitdistrと同じ結果が得られることが 分かり、floc制約をしない場合、どのような推定結果がでるのか調べて見ました。

今回、可視化はmatplotを使いました。floc=0の設定をしない方がより良くフィット していることが分かります。

# locも含めた分布のフィッティング結果 samp = np.array(davis.weight) fit_alpha, fit_loc, fit_beta = stats.gamma.fit(samp) print fit_alpha, fit_loc, fit_beta x = np.linspace(20,180,200) plt.figure() # fitted distribution pdf_fitted = stats.gamma.pdf(x, fit_alpha, loc=fit_loc, scale=fit_beta) plt.title('Gamma distribution') plt.plot(x,pdf_fitted,'r-') plt.hist(samp,bins=24, normed=1,alpha=.3) GgSave("5.png", 0.7) 
        
4.37074130051 36.4314864684 6.71934493135
Saving 8.0 x 6.0 in image.

4.37074130051 36.4314864684 6.71934493135
Saving 8.0 x 6.0 in image.
# fitdistrによるフィッティング結果 fit_alpha = 22.48548 fit_loc = 0 fit_beta = 2.926333 x = np.linspace(20,180,200) plt.figure() # fitted distribution pdf_fitted = stats.gamma.pdf(x, fit_alpha, loc=fit_loc, scale=fit_beta) plt.title('Gamma distribution') plt.plot(x,pdf_fitted,'r-') plt.hist(samp,bins=24, normed=1,alpha=.3) GgSave("6.png", 0.7) 
        
Saving 8.0 x 6.0 in image.

Saving 8.0 x 6.0 in image.

訓練データに異常データがまざっている場合

データに異常データが混ざっている場合の例として、以下の様に2種類の正規分布が混ざっている場合に ついて考えてみます。

# 訓練データに異常データがまざっている場合（2個の正規分布の混合分布を考える） x = var('x') N = 1000 mu0 = 3 mu1 = 0 sig0 = 0.5 sig1 = 3 pi0 = 0.6 pi1 = 0.4 # 分布を表示 g0 = plot(_gauss(x, mu0, sig0^2), [x, -5, 6], rgbcolor='green') g1 = plot(_gauss(x, mu1, sig1^2), [x, -5, 6], rgbcolor='blue') g01 = plot(_gauss(x, mu0, sig0^2)+ _gauss(x, mu1, sig1^2), [x, -5, 6], rgbcolor='red') (g0 + g1 + g01).show(figsize=5) 
       

分布は以下の様な混合正規分布として表されます。 $$ p(x) = \pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) + \pi_1 \mathcal{N}(x | \mu_1, \sigma_1^2) $$

パラメータを$ \theta = (\pi_0, \mu_0, \sigma_0^2, \pi_1, \mu_1, \sigma_1^2)$にまとめると、 対数尤度は、以下の様になります。 $$ L(\theta, \mathcal{D}) = \sum_{n=1}^N ln \left \{ \pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) + \pi_1 \mathcal{N}(x | \mu_1, \sigma_1^2) \right \} $$

これを$\mu_0$で偏微分すると、以下の様になります。 $$ 0 = \frac{\partial L}{\partial \mu_0} = \sum_{n=1}^N \frac{\pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) }{\pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) + \pi_1 \mathcal{N}(x | \mu_1, \sigma_1^2)} \frac{(x^{(n)} - \mu_0)}{\sigma_0^2} $$ 同様に$\sigma_0^2$で偏微分すると、以下の様になります。 $$ 0 = \frac{\partial L}{\partial \sigma_0^2} = \sum_{n=1}^N \frac{\pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) }{\pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) + \pi_1 \mathcal{N}(x | \mu_1, \sigma_1^2)} \frac{ \{-(x^{(n)} - \mu_0)^2 + \sigma_0^2 \}}{2} $$ ここで、データが$\pi_i$のどちらの集合に属するかその期待値を帰属度として以下の様に定義します。 $$ q_i^{(n)} = \frac{\pi_i \mathcal{N}(x | \mu_i, \sigma_i^2) }{\pi_0 \mathcal{N}(x | \mu_0, \sigma_0^2) + \pi_1 \mathcal{N}(x | \mu_1, \sigma_1^2)} $$

$\mu_i, \sigma_i^2$の推定値は、以下の様に求まります。 $$ \hat{\mu}_i = \frac{\sum_{n=1}^N q_i^{(n)} x^{(n)}}{\sum_{n'=1}^N q_i^{(n')}} $$ $$ \hat{\sigma}_i^2 = \frac{\sum_{n=1}^N q_i^{(n)} (x^{(n)} - \mu_i)^2}{\sum_{n'=1}^N q_i^{(n')}} $$

また、$\pi_i$は、以下の様になります。 $$ \hat{\pi}_i = \frac{1}{N} \sum_{n=1}^N q_i^{(n)} $$

EMアルゴリズム

混合正規分布をEMアルゴリズムを使って求めます。

EMアルゴリズムでは、以下の手順でパラメータの値を求めます。

  1. パラメータ$\theta$の初期値をセット
  2. 潜在変数（ここではどの正規分布に属するか）の期待値$q_i^{(n)}$が最大になるようなパラメータを求める
  3. 求まった潜在変数からパラメータの値を再設定する
  4. 値が収束するまで２）に戻る

テストデータの生成

以下の2つの正規分布をもつテストデータ1000個を生成します。 $$ (\pi_0, \pi_1) = (0.6, 0.4), (\mu_0, \sigma_0) = (3, 0.5), (\mu_1, \sigma_1) = (0, 3) $$

# テストデータの生成 T0 = RealDistribution('gaussian', sig0, seed=1) d0 = [T0.get_random_element()+mu0 for i in range(pi0*N)] h0 = histogram(d0, bins=20, normed=True, color='green') T1 = RealDistribution('gaussian', sig1, seed=1) d1 = [T1.get_random_element()+mu1 for i in range(pi1*N)] h1 = histogram(d1, bins=20, normed=True, color='blue') # プロットして確認 (h0+h1+g0+g1+g01).show(figsize=5) 
       

初期値のセット

各パラメータの初期値を以下の様に設定します。 $$ (\pi_0, \pi_1) = (0.5, 0.5), (\mu_0, \sigma_0) = (5, 1.0), (\mu_1, \sigma_1) = (-5, 5.0) $$

# データを結合 d = d0 + d1 # pi, sig, muの初期値 pi0 = 0.5 pi1 = 0.5 mu0 = 5.0 mu1 = -5.0 sig0 = 1.0 sig1 = 5.0 
       

EMアルゴリズムの実行

EMアルゴリズムを10回繰り返します。

高々10回の繰り返しで、とても精度良くパラメータが推定されています。

# EMアルゴリズムを10回繰り返す x = np.array(d) for i in range(10): piN0 = pi0*stats.norm.pdf(x, mu0, sig0); piN1 = pi1*stats.norm.pdf(x, mu1, sig1) qn0 = piN0/(piN0 + piN1); qn1 = piN1/(piN0 + piN1) pi0 = sum(qn0)/N; pi1 = sum(qn1)/N mu0 = sum(qn0*x)/(N*pi0); mu1 = sum(qn1*x)/(N*pi1) sig0 = sqrt(sum(qn0*(x - mu0)*(x - mu0))/(N*pi0)) sig1 = sqrt(sum(qn1*(x - mu1)*(x - mu1))/(N*pi1)) 
       
# 結果出力 print pi0, pi1 print mu0, mu1 print sig0, sig1 
        
0.602760455466 0.397239544534
2.98730105891 -0.292164057442
0.512590277559 2.85810740491
0.602760455466 0.397239544534
2.98730105891 -0.292164057442
0.512590277559 2.85810740491
# g0とg1を合わせて正規化したグラフ x = var('x') g01 = plot((_gauss(x, mu0, sig0^2)+ _gauss(x, mu1, sig1^2))/2, [x, -9, 9], rgbcolor='red') hd = histogram(d, bins=50, normed=True) (hd + g01).show(figsize=5) 
       
# 帰属度qn0, qn1の分布 q0 = list_plot(list(qn0), rgbcolor='red') q1 = list_plot(list(qn1)) (q0 + q1).show(figsize=5) 
       

パッケージを使って混合正規分布（GMM）を求める

最後にScikit-learnの混合正規分布パッケージsklearn.mixtureを使って、 davisの体重・身長のデータに対して、混合正規分布を求めてみます。

GMMは、異常なデータに影響を受けるため、12番目のデータを除いた集合で、 混合正規分布を求めます。

GMMの引数のn_componentsが求める分布の数です。RのMclustを使うとベストな分布数も計算できますが、 ここでは2として計算します。

# 混合正規分布 # SageでのGMMは、http://www15191ue.sakura.ne.jp:8000/home/pub/57/ を参照してください。 from sklearn import mixture # 分類器の生成 classifier = mixture.GMM(n_components=2, covariance_type='full') 
       
# GMMを求める # 12番目のデータを除く（体重と身長を間違えたデータ） davis_e11 = davis[davis.index != 11] data = davis_e11[['weight', 'height']].values # 混合正規分布を求める classifier.fit(data) 
        
GMM(covariance_type='full', init_params='wmc', min_covar=0.001,
  n_components=2, n_init=1, n_iter=100, params='wmc', random_state=None,
  thresh=None, tol=0.001)
GMM(covariance_type='full', init_params='wmc', min_covar=0.001,
  n_components=2, n_init=1, n_iter=100, params='wmc', random_state=None,
  thresh=None, tol=0.001)

分類結果のプロット

大切なのは、GMMによる分類ではなく、その結果を可視化して特徴を把握することです。

GMMによってデータがどのグループに属するか分類した結果は、predict関数で取得し、 davis_e11のpredにセットします。

次にggplotでpredで色を変えて分布をプロットします。 たったこれだけで、GMMの結果を確認することができます。

# 分類結果をpredにセット davis_e11['pred'] = classifier.predict(data) # 分類結果の表示 # 身長（height）と体重（weight）の関係をみる ggplot(davis_e11, aes(x='weight', y='height', colour='pred')) + geom_point() 
        
<ggplot: (8738328197473)>
<ggplot: (8738328197473)>
GgSave('fig_3.8.png', fac=0.8) 
        
Saving 11.0 x 8.0 in image.

Saving 11.0 x 8.0 in image.

正規分布の形もプロット

GMMの結果は、$\mu$がmeans_、$\sigma$がcovars_、$\pi$がweights_変数にセットされています。

Sageのcontour_plot関数を使って正規分布をプロットしてみます。

# クラス分けされたデータ cls0 = davis_e11[davis_e11.pred == 0][['weight', 'height']].values cls1 = davis_e11[davis_e11.pred == 1][['weight', 'height']].values 
       
# 正規分布の推定値 mu0 = vector(classifier.means_[0]) mu1 = vector(classifier.means_[1]) sig0 = matrix(classifier._get_covars()[0]) sig1 = matrix(classifier._get_covars()[1]) view(mu0); view(mu1) view(sig0); view(sig1) 
       
# ベクトルのガウス分布を定義 def _gauss(v, mu, sigma): d = len(v); sigma_inv = sigma.inverse(); sigma_abs_sqrt = sigma.det().sqrt(); val = -(v - mu) * sigma_inv * (v - mu).column()/2; a = (2*pi)**(-d/2) * sigma_abs_sqrt**-0.5 return a * e**val[0]; 
       
x, y = var('x y') # データのプロット pt_plt = Graphics() for pt in cls1: pt_plt += point(pt, rgbcolor='red') for pt in cls0: pt_plt += point(pt, rgbcolor='blue') # 正規分布を表示 r_cnt_plt = contour_plot(lambda x, y : _gauss(vector([x, y]), mu1, sig1), [x, 20, 130], [y, 140, 200], contours = 1, cmap=['red'], fill=False) b_cnt_plt = contour_plot(lambda x, y : _gauss(vector([x, y]), mu0, sig0), [x, 20, 130], [y, 140, 200], contours = 1, cmap=['blue'], fill=False) (pt_plt + r_cnt_plt + b_cnt_plt).show(figsize=5) 
       

混合正規分布データの異常度

混合正規分布の異常度は、以下の様になります。 $$ a(x') = -ln \left \{ \sum_{k=1}^K \hat{\pi}_k \mathcal{N}(x' | \hat{\mu}_k, \hat{\Sigma}_k) \right \} $$

#異常度の計算 # pi_nを求める pi_vec = vector(classifier.weights_); view(pi_vec) # davis全データに対する事後分布を求める a = [- log(pi_vec[0]*_gauss(vector(x), mu0, sig0)+pi_vec[0]*_gauss(vector(x), mu1, sig1)).n() for x in davis[['weight', 'height']].values ] 
       
list_plot(a, figsize=5)